{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**NOTE** This notebook is not part of the documentation. It's not meant to be in the webpage. It's something I wrote when I was testing the new functionality and I think it's nice to have it handy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import bambi as bmb\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = bmb.load_data(\"sleepstudy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Reaction</th>\n",
       "      <th>Days</th>\n",
       "      <th>Subject</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>249.5600</td>\n",
       "      <td>0</td>\n",
       "      <td>308</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>258.7047</td>\n",
       "      <td>1</td>\n",
       "      <td>308</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>250.8006</td>\n",
       "      <td>2</td>\n",
       "      <td>308</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>321.4398</td>\n",
       "      <td>3</td>\n",
       "      <td>308</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>356.8519</td>\n",
       "      <td>4</td>\n",
       "      <td>308</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Reaction  Days  Subject\n",
       "0  249.5600     0      308\n",
       "1  258.7047     1      308\n",
       "2  250.8006     2      308\n",
       "3  321.4398     3      308\n",
       "4  356.8519     4      308"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "       Formula: Reaction ~ 1 + Days + (1 + Days | Subject)\n",
       "        Family: gaussian\n",
       "          Link: mu = identity\n",
       "  Observations: 180\n",
       "        Priors: \n",
       "    target = mu\n",
       "        Common-level effects\n",
       "            Intercept ~ Normal(mu: 298.5079, sigma: 261.0092)\n",
       "            Days ~ Normal(mu: 0.0, sigma: 48.8915)\n",
       "        \n",
       "        Group-level effects\n",
       "            1|Subject ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 261.0092))\n",
       "            Days|Subject ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 48.8915))\n",
       "        \n",
       "        Auxiliary parameters\n",
       "            sigma ~ HalfStudentT(nu: 4.0, sigma: 56.1721)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = bmb.Model(\"Reaction ~ 1 + Days + (1 + Days | Subject)\", data)\n",
    "model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [sigma, Intercept, Days, 1|Subject_sigma, 1|Subject_offset, Days|Subject_sigma, Days|Subject_offset]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "862019c9e33b4cf18f6acc16c95cd34e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 5 seconds.\n"
     ]
    }
   ],
   "source": [
    "idata = model.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Reaction</th>\n",
       "      <th>Days</th>\n",
       "      <th>Subject</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>249.5600</td>\n",
       "      <td>0</td>\n",
       "      <td>xxx</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>258.7047</td>\n",
       "      <td>1</td>\n",
       "      <td>xxx</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>250.8006</td>\n",
       "      <td>2</td>\n",
       "      <td>xxx</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>321.4398</td>\n",
       "      <td>3</td>\n",
       "      <td>xxx</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>356.8519</td>\n",
       "      <td>4</td>\n",
       "      <td>xxx</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>414.6901</td>\n",
       "      <td>5</td>\n",
       "      <td>xxx</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>382.2038</td>\n",
       "      <td>6</td>\n",
       "      <td>xxx</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>290.1486</td>\n",
       "      <td>7</td>\n",
       "      <td>xxx</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>430.5853</td>\n",
       "      <td>8</td>\n",
       "      <td>xxx</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>466.3535</td>\n",
       "      <td>9</td>\n",
       "      <td>xxx</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>249.5600</td>\n",
       "      <td>0</td>\n",
       "      <td>308</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>258.7047</td>\n",
       "      <td>1</td>\n",
       "      <td>308</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>250.8006</td>\n",
       "      <td>2</td>\n",
       "      <td>308</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>321.4398</td>\n",
       "      <td>3</td>\n",
       "      <td>308</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>356.8519</td>\n",
       "      <td>4</td>\n",
       "      <td>308</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>414.6901</td>\n",
       "      <td>5</td>\n",
       "      <td>308</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>382.2038</td>\n",
       "      <td>6</td>\n",
       "      <td>308</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>290.1486</td>\n",
       "      <td>7</td>\n",
       "      <td>308</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>430.5853</td>\n",
       "      <td>8</td>\n",
       "      <td>308</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>466.3535</td>\n",
       "      <td>9</td>\n",
       "      <td>308</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    Reaction  Days Subject\n",
       "0   249.5600     0     xxx\n",
       "1   258.7047     1     xxx\n",
       "2   250.8006     2     xxx\n",
       "3   321.4398     3     xxx\n",
       "4   356.8519     4     xxx\n",
       "5   414.6901     5     xxx\n",
       "6   382.2038     6     xxx\n",
       "7   290.1486     7     xxx\n",
       "8   430.5853     8     xxx\n",
       "9   466.3535     9     xxx\n",
       "10  249.5600     0     308\n",
       "11  258.7047     1     308\n",
       "12  250.8006     2     308\n",
       "13  321.4398     3     308\n",
       "14  356.8519     4     308\n",
       "15  414.6901     5     308\n",
       "16  382.2038     6     308\n",
       "17  290.1486     7     308\n",
       "18  430.5853     8     308\n",
       "19  466.3535     9     308"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_new = data.head(10).reset_index(drop=True)\n",
    "df_new[\"Subject\"] = \"xxx\"\n",
    "df_new = pd.concat([df_new, data.head(10)])\n",
    "df_new = df_new.reset_index(drop=True)\n",
    "df_new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = model.predict(idata, data=df_new, inplace=False, sample_new_groups=True)\n",
    "\n",
    "reaction_draws = p.posterior[\"mu\"]\n",
    "mean = reaction_draws.mean((\"chain\", \"draw\")).to_numpy()\n",
    "bounds = reaction_draws.quantile((0.025, 0.975), (\"chain\", \"draw\")).to_numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAz8AAAF0CAYAAAANVYfFAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAbfVJREFUeJzt3Xl8VPW9P/7X7EtmMllnshAgQEAhoCxqQRSVxQW01bbYWu+tre23ty5Xrnj1Wnur9LZSvV/teqVf/Xm1V2rx9tZexYWSsEQQFAyLQAIhECAhmSXLbJn9nM/vj5CRkASyTzLzej4eGHLmzOQzE5z3vM5nUwghBIiIiIiIiJKcMtENICIiIiIiGgkMP0RERERElBIYfoiIiIiIKCUw/BARERERUUpg+CEiIiIiopTA8ENERERERCmB4YeIiIiIiFICww8REREREaUEhh8iIiIiIkoJDD80aJ988gm+/vWvIz8/H1qtFnl5efja176G3bt39+txnnnmGSgUigG1Yfv27VAoFNi+ffuA7t9XN9xwA2644YZh/Rn9FQgE8Mwzz/T43F9//XUoFAqcOnWqX485mN/FcNi1axeeeeYZuN3uRDeFiJJU5/ulXq/H6dOnu91+ww03oLS0NAEtSxzWF0pGDD80KL/97W9x7bXXoqGhAc8//zzKy8vxf//v/8XZs2excOFC/O53v+vzY33ve9/rd2DqNGfOHOzevRtz5swZ0P3HskAggDVr1vRYnJYvX47du3cjPz+/X485mN/FcNi1axfWrFnD4kREwy4cDuPHP/5xopsxKrC+UDJSJ7oBNHZ9/PHHWLVqFW677Tb89a9/hVr9xT+nb3zjG7jzzjvxyCOPYPbs2bj22mt7fZxAIACj0Yhx48Zh3LhxA2pLeno6vvSlLw3ovmOVEAKhUOii5+Tm5iI3N7ffjz2Y3wUR0Vh2yy234M0338Rjjz2GK664ItHNSQjWF0pm7PmhAVu7di0UCgXWrVvXJfgAgFqtxksvvQSFQoFf/OIX8eOd3d379u3D1772NWRmZmLy5MldbjtfOBzG6tWrkZeXB6PRiOuvvx6VlZWYOHEi7rvvvvh5PQ17u++++2AymVBbW4vbbrsNJpMJRUVFWL16NcLhcJefs2bNGlxzzTXIyspCeno65syZg1dffRVCiAG9NhMnTsSKFSvw17/+FbNmzYJer8ekSZPwm9/8pst5oVAIq1evxpVXXgmLxYKsrCzMnz8f77zzTrfHVCgUeOihh/D73/8el19+OXQ6Hf7whz/Ei8+aNWugUCigUCjir01vwxI2bdqExYsXw2KxwGg04vLLL8fatWvjt/f0uxjO5/TGG2/g8ssvh9FoxBVXXIH33nuvS1v++Z//GQBQXFwcf47DPcSRiFLT448/juzsbDzxxBOXPFcIgZdeeglXXnklDAYDMjMz8bWvfQ0nT56Mn/Mf//EfUCqVcDqd8WMvvPACFAoFHnzwwfgxWZaRmZmJ1atXX/Rnsr6wvtDgsOeHBkSSJGzbtg3z5s3r9QpOUVER5s6di61bt0KSJKhUqvhtd911F77xjW/gH/7hH9De3t7rz/nOd76Dt956C48//jhuuukmVFVV4c4774TX6+1TO6PRKO644w7cf//9WL16NT766CP827/9GywWC37yk5/Ezzt16hR+8IMfYPz48QA65jE9/PDDOHv2bJfz+uPAgQNYtWoVnnnmGeTl5eGPf/wjHnnkEUQiETz22GMAOsJda2srHnvsMRQWFiISiaC8vBx33XUXXnvtNfz93/99l8f83//9X+zYsQM/+clPkJeXh6ysLGzatAm33HIL7r//fnzve98DgItejXv11Vfx/e9/H4sWLcLvf/97WK1W1NTU4PDhwwl5Tu+//z727t2Ln/70pzCZTHj++edx55134tixY5g0aRK+973vobW1Fb/97W/x9ttvx4dYTJ8+ve+/DCKiPjKbzfjxj3+MRx55BFu3bsVNN93U67k/+MEP8Prrr+Mf//Ef8dxzz6G1tRU//elPsWDBAhw8eBA2mw1LliyBEAJbtmzBN7/5TQBAeXk5DAYDysrK4o/12Wefwe12Y8mSJZdsI+sL6wsNgiAaALvdLgCIb3zjGxc97+677xYAhMPhEEII8fTTTwsA4ic/+Um3cztv63TkyBEBQDzxxBNdzvvTn/4kAIhvf/vb8WPbtm0TAMS2bdvix7797W8LAOK///u/u9z/tttuE9OmTeu1zZIkiWg0Kn7605+K7OxsIcty/LZFixaJRYsWXfQ5CyHEhAkThEKhEAcOHOhyfOnSpSI9PV20t7f3eL9YLCai0ai4//77xezZs7vcBkBYLBbR2tra5bjL5RIAxNNPP93t8V577TUBQNTV1QkhhPD5fCI9PV0sXLiwy/O60IW/i+F8TjabTXi93vgxu90ulEqlWLt2bfzYv//7v3d5HkREQ63z/XLv3r0iHA6LSZMmiXnz5sXfKxctWiRmzJgRP3/37t0CgHjhhRe6PE59fb0wGAzi8ccfjx8bN26c+O53vyuEECIcDou0tDTxxBNPCADi9OnTQgghfv7znwuNRiP8fv9F28n60vfnxPpCPeGwNxpW4tywsQu7uL/61a9e8r4VFRUAgJUrV3Y5/rWvfa3bMLveKBQK3H777V2OzZo1q9tKPlu3bsWSJUtgsVigUqmg0Wjwk5/8BC0tLV2GKvTHjBkzuo0Xv+eee+D1erFv3774sT//+c+49tprYTKZoFarodFo8Oqrr6K6urrbY950003IzMwcUHuAjomdXq8XDzzwwIBW2xmO53TjjTfCbDbHv7fZbLBarT2utkRENBK0Wi1+9rOf4bPPPsN///d/93jOe++9B4VCgXvvvRexWCz+Jy8vD1dccUWXoVOLFy9GeXk5gI734UAggEcffRQ5OTnx3p/y8nLMnz8faWlpl2wf6wvrCw0cww8NSE5ODoxGI+rq6i563qlTp2A0GpGVldXleF9Wh2lpaQHQ8WZ1PrVajezs7D6102g0Qq/Xdzmm0+m6TOTcs2cPli1bBgB45ZVX8PHHH2Pv3r146qmnAADBYLBPP+tCeXl5vR7rfG5vv/02Vq5cicLCQqxfvx67d+/G3r178d3vfrfHyab9XVXnQi6XCwAGPNl0OJ5TT79LnU434NediGgofOMb38CcOXPw1FNPIRqNdrvd4XBACAGbzQaNRtPlzyeffILm5ub4uUuWLMGZM2dw/PhxlJeXY/bs2bBarbjppptQXl6OYDCIXbt29WnIG8D6wvpCg8E5PzQgKpUKN954IzZt2oSGhoYe3+waGhpQWVmJW2+9tct8H6B7T1BPOt+0HA4HCgsL48djsVj8jXAobNiwARqNBu+9916XoPS///u/g3pcu93e67HO57Z+/XoUFxfjrbfe6vKaXLggQ6fB7o3QOVa7oaFhQPcfjudERDQaKRQKPPfcc1i6dClefvnlbrfn5ORAoVBgx44d0Ol03W4//9jixYsBdPTulJWVYenSpfHjP/7xj/HRRx8hHA73OfywvrC+0MCx54cG7Mknn4QQAg888AAkSepymyRJ+OEPfwghBJ588skBPf71118PAHjrrbe6HP+f//kfxGKxgTW6BwqFAmq1uktACwaDeOONNwb1uEeOHMHBgwe7HHvzzTdhNpvj+xEpFApotdoub+J2u73HlWt601lg+3Ila8GCBbBYLPj9738/oJXsRuo5Xag/z5GIaKgsWbIES5cuxU9/+lP4/f4ut61YsQJCCJw9exbz5s3r9mfmzJnxc/Pz8zF9+nT85S9/QWVlZTz8LF26FC6XCy+++CLS09Nx1VVX9aldrC+sLzRwDD80YNdeey1+9atf4f3338fChQvxxz/+ETt27MAf//hHXHfddfjggw/wq1/9CgsWLBjQ48+YMQPf/OY38cILL+BHP/oRysvL8etf/xqPP/44LBYLlMqh+ee7fPly+P1+3HPPPSgrK8OGDRtw3XXX9Xglrz8KCgpwxx134LXXXsOmTZtw7733oqysDD/+8Y9hNBoBdBTPY8eO4YEHHsDWrVvxhz/8AQsXLuzX8AOz2YwJEybgnXfewebNm/HZZ5/1uuO2yWTCCy+8gI8++ghLlizBhg0bsG3bNrzyyit46KGHRs1zulDnh4hf//rX2L17Nz777DP4fL4BPx4RUV8999xzcLlcqKys7HL82muvxf/5P/8H3/nOd/D444/jvffew7Zt2/Dmm2/igQcewLp167qcv3jxYmzZsgVarTa+911xcTGKi4uxefNm3HDDDX2ez8r6wvpCg5DAxRYoSezevVt87WtfEzabTajVamG1WsVdd90ldu3a1e3czlVeXC5Xr7edLxQKiUcffVRYrVah1+vFl770JbF7925hsVjEP/3TP8XP6221t7S0tD79nP/8z/8U06ZNEzqdTkyaNEmsXbtWvPrqq91WgOnPam/Lly8X//M//yNmzJghtFqtmDhxonjxxRe7nfuLX/xCTJw4Ueh0OnH55ZeLV155pcc2AhAPPvhgjz+vvLxczJ49W+h0ui4r4V24Gk+nDz74QCxatEikpaUJo9Eopk+fLp577rmLvkYj+ZwmTJjQZTU/IYR48sknRUFBgVAqld1+10REg3X+am8XuueeewSALqu9dfrP//xPcc0114i0tDRhMBjE5MmTxd///d+Lzz77rMt577zzjgAgli5d2uX497//fQFA/OY3v+lTO1lfWF9ocBRCDHAXR6IE2bVrF6699lr88Y9/xD333JPo5vRo4sSJKC0t7bKZ2liXjM+JiGisScb34mR8TjR6ccEDGtXKysqwe/duzJ07FwaDAQcPHsQvfvELlJSU4K677kp084iIiIhoDGH4oVEtPT0dmzdvxq9+9Sv4fD7k5OTg1ltvxdq1a7stYU1EREREdDEc9kZERERERCmBq70REREREVFKYPghIiIiIqKUwPBDREREREQpYUwueCDLMhobG2E2m7vs8ktERMNPCAGfz4eCgoIh22w4GbA2ERElRn/q0pgMP42NjSgqKkp0M4iIUlp9fT3GjRuX6GaMGqxNRESJ1Ze6NCbDj9lsBtDxBNPT0xPcGiKi1OL1elFUVBR/L6YOrE1ERInRn7o0JsNP53CC9PR0FhgiogTh0K6uWJuIiBKrL3WJg7WJiIiIiCglMPwQEREREVFKYPghIiIiIqKUwPBDREREREQpgeGHiIiIiIhSAsMPERERERGlhDG51DUREQ2MJAvsqWuF0xeC1azH1cVZUCm5ZDURESXGSNclhh8iohSx6XAT1mysQpMnFD+Wb9Hj6dun45bS/AS2jIiIUlEi6hKHvRERpYBNh5vww/X7uhQYALB7Qvjh+n3YdLgpQS0jIqJUlKi6xPBDRJTkJFlgzcYqiB5u6zy2ZmMVJLmnM4iIiIZWIusSww8RUZLbU9fa7cra+QSAJk8Ie+paR65RRESUshJZlxh+iIiSnNPXe4EZyHlERESDkci6xPBDRJTkrGb9kJ5HREQ0GImsSww/RERJ7uriLORb9Oht4VAFOlbXubo4aySbRUREKSqRdYnhh4goyamUCjx9+3QA6FZoOr9/+vbp3O+HiIhGRCLrEsMPEVEKuKU0H+vunYM8S9chBHkWPdbdO4f7/BAR0YhKVF3iJqdERCniltJ8LJ2eN6I7aRMREfUmEXWJ4YeIKIWolArMn5yd6GYQEREBGPm6xGFvRERERESUEhh+iIhSTF1zO+qa2xPdDCIiorhgREJMkof953DYGxFRigjHJGw/6kR5tRPTC9JRnFOc6CYREVGKC0YkfFLXgp3Hm3FraR7mTRzebRcYfoiIUkB9awDvHmjEobMehGISZCES3SQiIkphMUnGwQY3yqsdqHMFEJVkhKLs+SEiokGISTI+PtGCvx22ozUQwcRsIxo9oUQ3i4iIUpQQAsedfpRXOVDV5IVKqcCk3DScGqHh2Aw/RERJyukN4d2Djdh3ug1pOjWmWk1QKLisNRERJUaTJ4gt1U5Unm5DOCqhMNMAo3Zk48igFjxYu3YtFAoFVq1aFT923333QaFQdPnzpS99qcv9wuEwHn74YeTk5CAtLQ133HEHGhoaBtMUIiI6R5YF9tS14qXttdh7qhWFmQYUZBgYfIiIKCE8wSjeO9iIX5cfx44aFywGDUps5hEPPsAgen727t2Ll19+GbNmzep22y233ILXXnst/r1Wq+1y+6pVq7Bx40Zs2LAB2dnZWL16NVasWIHKykqoVKqBNomIKOW5AxG893kTPjnZAo1Kiak2M5QMPURElAChqITPTrVhy1EHGt1BZKVpMS3PnNCLcQMKP36/H9/61rfwyiuv4Gc/+1m323U6HfLy8nq8r8fjwauvvoo33ngDS5YsAQCsX78eRUVFKC8vx8033zyQJhERpTQhBA6d9WDjwUacbglgXKYBZr0m0c0iIqIUJMsdNams2oETTj8MGhVKrGaolIm/GDegYW8PPvggli9fHg8vF9q+fTusViumTp2K73//+3A6nfHbKisrEY1GsWzZsvixgoIClJaWYteuXT0+Xjgchtfr7fKHiIg6+MMx/GVfA17dWQeHN4SpNjODzwhgbSIi6koIgZMuP/6/nSfx/+04iVPN7ZiQbURRlnFUBB9gAD0/GzZswL59+7B3794eb7/11lvx9a9/HRMmTEBdXR3+9V//FTfddBMqKyuh0+lgt9uh1WqRmZnZ5X42mw12u73Hx1y7di3WrFnT36YSESW9GocP7xw4i+NOP/LS9cg0ai99JxoSrE1ERF9w+kLYfsyFT0+2IBCRUJhhQJpu9K2t1q8W1dfX45FHHsHmzZuh1+t7POfuu++O/720tBTz5s3DhAkT8P777+Ouu+7q9bGFEL2O/3vyySfx6KOPxr/3er0oKirqT9OJiJJKKCphS7UDW486EY7JmJJrgkY1qDVsqJ9Ym4iIOkYffFzbjIpjLjT7w8hL16NwFC+y06/wU1lZCafTiblz58aPSZKEjz76CL/73e8QDoe7LViQn5+PCRMm4Pjx4wCAvLw8RCIRtLW1den9cTqdWLBgQY8/V6fTQafT9aepRERJ63RLO9450IgjjR7kmHSjusgkM9YmIkplkZiMytNt2HrUgTOtAWQYOhYzGO2L7PQr/CxevBiHDh3qcuw73/kOLrvsMjzxxBM9rtTW0tKC+vp65OfnAwDmzp0LjUaDsrIyrFy5EgDQ1NSEw4cP4/nnnx/o8yAiSnpRScaOGhf+VuWANxhFcU4adGqukElERCNHlgWqmrzYUu3AMbsPWrUSU3JNUI+R0Qf9Cj9msxmlpaVdjqWlpSE7OxulpaXw+/145pln8NWvfhX5+fk4deoUfvSjHyEnJwd33nknAMBiseD+++/H6tWrkZ2djaysLDz22GOYOXNmrwsoEBGluiZPEO8caMTBejfS9RqUcMNSIiIaYWdaAthy1IH9Z9ogyUBRlhF6zdi6CDeks5BUKhUOHTqE//qv/4Lb7UZ+fj5uvPFGvPXWWzCbzfHzfvnLX0KtVmPlypUIBoNYvHgxXn/9de7xQ0R0AUkW+LSuBR8eaoLTF8aErDQYtHyvJCKikdPaHsH2Y07sPtECbyiKwoyxu53CoMPP9u3b4383GAz429/+dsn76PV6/Pa3v8Vvf/vbwf54IqKk1eIPY+Pnjdhb1wa9RolptsRuDEdERKklEInhk5Mt2HbUBYc3BKtZN+Zr0ehbf46IKMUJIbC/3o33DjaioS2IoiwjTKNwuVAiIkpOMUnGgXo3yqsdONXcDrNeg2k2M5SjZK+ewWA1JSIaRbyhKD74vAm7TjRDAQWm2kbHjthERJT8hBCocfhRXuVAVZMXKqUCk5JsKwWGHyKiUUAIgeomH949eBYnXe0oyDDAYhib46mJiGjsaXQHsaXagcrTbYjEZIzLNCblHFOGHyKiBAtEYiircqDimAsRSUaJdewsGUpERGObJxDFjuMu7DjugjsYRb4luS++MfwQESXQSZcf7xxoxFG7F7kmHYpMxkQ3iYiIUkAoKmFPXSu2HnWi0R1EjmnsL2bQFww/REQJEI5JqDjmQnm1A/5QDJNyTNCq2dtDRETDS5YFDp31oKzagVqnH2laVUrNL2X4ISIaYQ1tAbx7oBGfN3iQYdRgCjcsJSKiEXCquR1lVQ583uCGAFCcnZZyF94YfoiIRkhMkrH7ZAs+PGRHa3sEE7LH3s7YREQ09rT4w9h2zIlPTrTAH5EwLsOAtBTdQiE1nzUR0Qhz+kJ472AjPjvthlGrwlQbe3uIiGh4BSMSdp9sxrajTji8YdjS9SjIMKR0/WH4ISIaRtGYjPWfnMa2Y04EIxJmjbPApE/eVXSIiCjxJFngQH0byqsdOOlqR/oFm5QKIXDWHUR7OIY0nRqFKRSIGH6IiIbJ/3xWj5++VwVvKBY/Vm33YtHUXEyxmhPYMiIiSkZCCJxwtaOsyo7DZ71QqxSYfMEmpbVOHypqXPCHpfgxk06VMrWJ4YeIaIjJssDvK2rx/N9qut3mD0t4/5Ady2ciJYoMERGNDKcvhK3VTuypa0UwKmFcpgFGbdeP+rVOH94/ZO9231SqTQw/RERDyB2IYOPnjfjNltqLnldR48LkXM77ISKiwfGHY/i4thkVx1xo9oeRZ9FjXGb3YWxCCFTUuC76WKlQmxh+iIiGgBACnzd48N7njdh/xo1QTL7o+f6whLPuIMZlclNTIiLqv6gkY9/pjnk9Z1oDyDBoMS3PDGUvweWsO9hlqFtPUqE2MfwQEQ2SJxjF347Y8fHxZggIZBn7tqBBezh26ZOIiIjOI4RAjcOPzVV2VDd5oVUpMSXXBLXq4vv19LXmJHttYvghIhogIQSONHqx8fNG1LnakW/RI8OoRUNboE/3T9U9FoiIaGCaPEFsqXbis1OtiEgyxmUaYejjfnF9rTnJXpuS+9kREQ0TfziGvx1uwo7jzZBkgRLrF1fdCjMMMOlUFx1eYNKpUJhhGKnmEhHRGOYNRbGjxoUdx5vR1h5BfoYBFkP/tk1gberA8ENE1A9CCBy1+7DxYCNqXX7kpeuRadR2OUehUGDR1NweV9TptGhqblJPKCUiosGLxGR8dqoV5dUOnHUHkWXsmNczkPrB2tSB4YeIqI/awzGUVztQccyFSEzGlAv2TjjfFKsZy2cipfdSICKigekcVr2l2oFjdh90GhWmWE1QKy8+r+dSWJsYfoiI+uS4w4d3DjaixuGD1aRDUdalV8KZYjVjcq4pZXfRJiKi/mtoC6C8yoH9Z9yIyQLjs4zQ9XFeT1+kem1i+CEiuohgRMKWage2HXMiFL14b09PFApFUi8ZSkREQ8MdiKCixoWPa5vhCUZRkGFAur5/83r6KpVrE8MPEVEvTrj8ePdAI442eZFt0qXUlTEiIhoZoaiEPXWt2HrUiUZ3EDkmHabZBjavhy6N4YeI6AKhqISKY06UVzvRHo5hUq4JWvXgxlkTERGdT5YFDp31oKzagVqnH2laFabazFApGXqGE8MPEdF5Tre0450Djahq9CAzTYspVhOvvhER0ZA61dyO8moHDta7IQBMzDZCpx66eT3UO4YfIiIA4ZiEHTUulFU54A3FMDEnjYWIiIiGVGt7BNuOOrD7ZCv8oRgKMw0wJfmmoqMNX22iMUiSBfbUtcLpC8Fq1uPq4ix2kw9CfWsA7x5sxKGzHlgMGpSwt4eIqN9Ym3oXjEj45GQLth51wuENwWrWocDGWpMIDD9EY8ymw01Ys7EKTZ5Q/Fi+RY+nb5+OW0rzE9iysScqyfi4thl/O2xHWzCKCVlG6IdwOVEiolTB2tSzmCTjYIMb5dUO1LnaYdZrMM1mhpKhMGE4g5doDNl0uAk/XL+vS3EBALsnhB+u34dNh5sS1LKxp9EdxKs76/DW3npIQmCq1cTgQ0Q0AKxN3QkhUOv04ZUdJ/GfO0+h0R3CpFwTCjIMDD4Jxp4fojFCkgXWbKyC6OE2AUABYM3GKiydnsdhBhcRk2TsPtmCTYftaPaHMSErDQYtQw8R0UCwNnXn8IawtdqJvadaEYxKGJdpgFHLj9yjBX8TRGPEnrrWblfVzicANHlC2FPXivmTs0euYWOIwxvCuwcasf9MG4w6NfdRICIaJNamL/hCUeyqbUZFTTOa/WHkW/QoykrNjURHM4YfojHC6eu9uAzkvFQiyQKf1rXgw0N2OH0hFGUakcbVdYiIBo21CYjEZFSebsOWagfOtAWQZdTisjxeXButWP2JxgirWT+k56UKly+M9z5vxGen2qDXKDHVZoaSBYmIaEikcm0SQqCqyYvyKgeO2X3QqpUoyTVBreKU+tGM4YdojLi6OAv5Fj3snlCPY6sVAPIsHUuLUsfO2Z+dbsMHh5rQ5AliXKaReykQEQ2xVK1NDW0BlFc7sP+MGzFJoIirhY4ZjKZEY4RKqcDTt08H0FFMztf5/dO3T0+ZCaUX0+wPY/2np/Ffu0/BE4xiqtXM4ENENAxSrTa5AxG8c+AsfrPlOHafaEVWmhZTuFromMLwQzSG3FKaj3X3zkGepevwgTyLHuvunZPSeykA5+b2nGzB77bWYufxZljNOozPMnJZUSKiYZQKtSkUlfBRjQu/Kj+O9z9vglqpxDSbCel6TaKbRv3ES6FEY8wtpflYOj2Pu2hfwOUL4/1DHXN7NColN5EjIhpByVqbJFng0FkPyqscqHX6YdKpMdVmHvPPK5Ux/BCNQSqlIumXDO0rSRbYU9eKvx1pQqO7YyU3k55vbUREIy2ZapMQAnXN7SivduDzBg8UAIpz0qBVc9DUWMdPCEQ0Zjm9Ibx/qAmVp9ugZW8PERENAZcvjO3HnPjkZAsCEQmFGQZuj5BE+JskojGnc9+eTYftcHhCKMrivj1ERDQ47eEYdp1oxvZjLjT7wrCl61GYYeB+PUmGnxaIaExxeEN47/NG7DvdBr1Ghal53LeHiIgGLibJ2F/vRnmVA6daArAY1KwtSYzhh4jGhJgk49O6Vmw63ASnL4yiTPb2EBHRwAkhUOPwo7zKgaomD9QqJabkpnGT0iQ3qN/u2rVroVAosGrVqvgxIQSeeeYZFBQUwGAw4IYbbsCRI0e63C8cDuPhhx9GTk4O0tLScMcdd6ChoWEwTSGiJNbkCeK1j0/hj5+eQTAiY6rNzOBDREQD1uQJYv0np7Fuey2ONHpQmGnExGwGn1Qw4N/w3r178fLLL2PWrFldjj///PN48cUX8bvf/Q579+5FXl4eli5dCp/PFz9n1apV+Otf/4oNGzZg586d8Pv9WLFiBSRJGvgzIaKkE5Vk7Djuwu+21uKz020ozNCjMNPAoQiDEIxIaI/EoOy2HSERUfLzBKN4//NG/Lr8OHbWNiPDqEWJzQwDNylNGQMKP36/H9/61rfwyiuvIDMzM35cCIFf/epXeOqpp3DXXXehtLQUf/jDHxAIBPDmm28CADweD1599VW88MILWLJkCWbPno3169fj0KFDKC8vH5pnRURjXqM7iFd31uFPn55BOCZjms0Eo5a9PQPlD8dw3OlHoyeIWYUWzJ+SHMvREhH1RTgm4ePaZvy6vAbvHGyEQgFMs5lhMXCT0lQzoPDz4IMPYvny5ViyZEmX43V1dbDb7Vi2bFn8mE6nw6JFi7Br1y4AQGVlJaLRaJdzCgoKUFpaGj+HiFJXVJJRUdPR27P/TBsKMw1cbWeAhBDwBqOocfjQ7AtjZmE6/mHRZDx44xRclpee6OYREQ07WRY4fNaDl7adwPpPTsMdiKLEaoLVrGddSVH9voy6YcMG7Nu3D3v37u12m91uBwDYbLYux202G06fPh0/R6vVdukx6jyn8/4XCofDCIfD8e+9Xm9/m01EY0BDWwAbDzbi8wYPTDo1ptnMLE4DIISAOxiFwxuCUavGVRMzce2UXJRYTdwHaQixNhGNbvWtAZRVO3DgTBtkAUzINkKn5vC2VNev8FNfX49HHnkEmzdvhl6v7/W8Cz+sCCEu+QHmYuesXbsWa9as6U9TiWgMicRkfFzbjM1H7GgLRDE+28jx1wMghEBLewTNvjDMBg2uK8nFgsnZKM5JY4gcBqxNRKOTJxBFRY0LO2td8ASjKMwwwKzn8Dbq0K9hb5WVlXA6nZg7dy7UajXUajUqKirwm9/8Bmq1Ot7jc2EPjtPpjN+Wl5eHSCSCtra2Xs+50JNPPgmPxxP/U19f359mE9EoVt8awP+34yT++7N6SEJgqs3E4NNPsizg8IZw1O5DTJKxdIYNq5aU4FvXjMekXBODzzBhbSIaXSIxGbtqm/GrLTV47/NGqJVKTLOZGXyoi371/CxevBiHDh3qcuw73/kOLrvsMjzxxBOYNGkS8vLyUFZWhtmzZwMAIpEIKioq8NxzzwEA5s6dC41Gg7KyMqxcuRIA0NTUhMOHD+P555/v8efqdDrodLp+PzkiGr3CMQk7alwor3bCHYxiQpYReoaefpFkAacvBHcgilyzDitmFeCaSVmwpffeM09Dh7WJaHQQQqC6yYeyKjuO2X0dG2DbzFBxmC/1oF/hx2w2o7S0tMuxtLQ0ZGdnx4+vWrUKzz77LEpKSlBSUoJnn30WRqMR99xzDwDAYrHg/vvvx+rVq5GdnY2srCw89thjmDlzZrcFFIgoOZ1uacfGg0040uhBul6DqVb2TvRHTJJh94bgD8dgS9fjq3OsuKo4G1lp2kQ3jYhoRDV5giivcuCz022ISQLjs4zQ8UIaXcSQrxv7+OOPIxgM4oEHHkBbWxuuueYabN68GWazOX7OL3/5S6jVaqxcuRLBYBCLFy/G66+/DpWK/1iJklkoKuGjGhe2VDvgDcUwnr09/RKJyWjyBBGKyijI0GPFFQWYOyET6RzSQUQpxh+OYUeNCxU1LrQFIiiwGJDOZaupDxRCCJHoRvSX1+uFxWKBx+NBejqXayUaC+qa2/HewUYcafQiw6iB1axjb08fhaISGj1BSJLA+GwjrivJxezxGQnb94jvwT3j60I0/GKSjH1n3CirsuNMawCZRi3rSZI47vDhnmsmYGFJTr/v25/3X+4YSETDKhSVUHHMiS1HnfCFYpiYzSEJfdUejsHuCUEAmJSbhuun5mJmoYW9ZUSUcoQQqHX6sfmIA0eaPNCqlJiSa4JaNaAtKymFMfwQ0bA54fLjvYONqG7yItOoRQnn9vSJLxSF3ROCSqnAtHwzrivJxYyCdGhY5IkoBTl9IWytdmJPXStCUQnjMo0waHkRiAaG4YeIhlwoKmFrtRPbjjnRHo5hYk4aN5a7BCEEvKEYHN4QdGolrijKwMKSHFyWl84Vi4goJQUiMeyqbca2Yy64fGHkW/QoyjImulk0xjH8ENGQqnX6sfFgI47avchO02EKe3suSgiBtkAUTl8YaToVvjQpG9dOycZk7s9DRClKkgUO1LtRXm3HSVc7LAYNpuWZoeR7Ig0Bhh8iGhKBSAxbq53YXuNCIBzDpBwTtGoO0+qNLARa/BG0+MOwGDW4cVou5k/OxvgsI0MPEaWsuuZ2lFXZ8XmDByqlApNzTRzyS0OK4YeIBkUIgSONXnx4uAm1znZkp2lRwN6eXkmygMsfRlt7BFlpWtxcmocvTcpGQYYh0U0jIkqY1vYIth11YPeJFrRHJBRmGJCm48dUGnr8V0VEA9baHsHmI3Z8crIFMVlgcm4ar9D1IibLcHjD8IWiyDXp8OUrC3B1cTZyzbpEN42IKGFCUQmf1rViS7UDdk8ItnQ9CjIMvIBGw4bhh4j6LSbJ2HuqDZuP2HHWHUReuh7tkRhOuvxI06lRyMIVF5Vk2D0hBCIx5FkMuHmGDVdNzEKGUZvophERJYwsCxxu9GBzlQO1Dj9MOjWm2cxQDuECL0IInHUH0R6OsTZRHMMPEfVLQ1sAHxxqwsF6N7RqFVRK4O39DfCHpfg5Jp0Ki6bmYorVnMCWJlYkJqPJE0QoJqMo04AvlxRizvgMmPXcgZyIUlt9awBl1Q7sP90GAaA4J23I54jWOn2oqHGxNlE3DD9E1CehqISPalzYetSJtkAE47OMaHQH8eFhR7dz/WEJ7x+yY/lMpFyRCUUlNHmCiEkCE3KMuG5KLmaPz+SeFESU8jyBKCpqXNhZ64InGMW4DCNM+qH/KFrr9OH9Q/Zux1O5NtEXGH6I6KKEEDhq9+HDQ02ocfiQYdRimq2jaFTUuC5634oaV8os2RyIxNDkCQEAJuWk4bqpuZg1zsL9jYgo5UViMj471Yryagca2oLIMekwzWYeltoghGBtooti+CGiXrkDEZRVObDrRAsiMRnF5y1f3dAW6DKcoCf+sISz7iDGZSbvpnT+cAx2TwgKBTDVZsZ1JTkoLbRw4QciSnlCCFQ3+VBWZccxuw96jQpTbeZh3bj5rDvI2kQXxfBDRN1IskDl6Tb87UgT6luDyLPokXnBBP32cKxPj9XX88YabzAKuzcErVqJmeMsWDglB5fnpw9rUSciGiuaPEGUVznw2ek2SLLA+CwjdJrh7wlP9dpEl8bwQ0RddMzjacL+M26olYper9L1df+FZNqnQQgBdzAKpzcMg1aFqyZm4topuSixmoZ0hSIiorHKH45h53EXth9zoS0QQYHFgHTDyC30koq1ifqHv3kiAtAxUX/ncRe2HnWhpT2CosyLbzBXmGGASae66PACk06FwiTYvFMIgdb2CFz+MEw6Na6dko0FU3IwKSeNY8aJiNCxrP++020or3bgTGsAmefmh470e2Qq1SYaGIYfohQnhMBxpx/vf96EY3YfLAYNptkuPRFUoVBg0dTcHlfU6bRoau6YDgeyEGjxR9DsDyPDqMGN06y4dkoOxmVyrwgiIqCjhtQ6/dh8xIEjTR5oVUpMyTVBnaB5j6lQm2hwGH6IUpgnGMWWagd21jYjFJUwMcfYr9XJpljNWD4TSbeXgiwLOP1htLVHkJWmxa2lefjS5GzkW3ilkIiok9MXwtZqJ/bUtSIUlTAu0zgqlvVP1tpEQ4PhhygFybLA/vo2bDpsx+mWAGzpehQNcNWbKVYzJueakmIXbUkWcHhD8ISisJp0uOOKAlwzKRu5Zl2im0ZENGoEIjF8XNuM7cdccPnCyLfoUZQ1ulZOS6baREOL4YcoxTi8IXx4qAmfnW6DSqFAic0EtXJwwxMUCsWYXjI0Jsmwe0PwhyXkpeuwZLoNV0/MQmaa9tJ3JiJKEZIscKDejfJqO0662pFh0GJanhnKURooxnptouHB8EOUIsIxCbtPtKCsyoFmf3jYdtYeSyIxGU2eIEJRGQUZeqy4ogBzJ2QiXT9yKxMREY0Fdc3t+NsROw6f9UClVGByron7mdGYlNqffIhSRK3Tjw8ONaGqyQuzTp2QFXhGk3BUQqMnhKgkY3y2EddNycWcCRkwavmWSER0vhZ/GNuOObH7RAsCEQmFGRdfCZRotOO/XqIk5gtFsaXaiZ3HXQhEJEwcoU3mRqtgREKjJwghBCbmpOG6klxcWZQBfQq/JkREPQlFJXxysgVbjzph94RgS9dzzgwlBYYfoiQkywIHG9zYdMSOuuZ2WE06FKRw0WoPx9DkCUGhAKZYTbiuJBczCy3Qqjlkg4jofLIscOisB2XVDtQ6/DDpO0YLcCNnShYMP0RJxukLYdNhO/bWtUKhAEoSuN9CovlCUdg9IahVCswoSMfCkhxMz09P2deDiOhi6lsDKKuyY/8ZNwCgOCeNF4ko6TD8ECWJSEzGJydbUFZlh9MXRmGGAeYUnbjvDUZh94agUysxe3wmFpbk8MolEVEvPIEottc48XFtMzzBKMZlGmHivB5KUvyXTZQE6prb8cHnTTjc6EGaVo2pttG79Ohw6gw9eo0SVxdn4bqSHEzONaXscD8ioosJxyTsrWtDebUDje4gcs26lF8Qh5Ifww/RGNYejmHbUee5XaxjGJ9lTMnJ+95zw9v0GiWumpiJ66fmMvQQEfVCCIGqJi82H3GgxuGDUavCVJsZKvaOUwpg+CEag2RZ4POzHvztSBNOuNqRk6ZDiTX1Pux3hh6dmqGHiKgvGt1BlFU5sO9MG2KSwIRsI3Tq1LtoRqmL4YdojGnyBPG3w3ZUnm7rWL0sBTeaY+ghIuofbyiKHTUufHS8Ge5ABAUWA9INqTkvlFIbww/RGBGMSNhZ68K2oy60tKfmgga+UBRN54We60pyMSUFe7yIiPoqKsmoPN2G8ioH6lsDyDJpOa+HUhrDD9EoJ4TAkUYvNh2247jThwxD6hWu80PP3AmZWDSVoYeI6GKEEKhx+LG5yo7qJh+0KgWm2ExQK1NrpADRhRh+iEYxpy+EzUcc2HuqFZIsMCnHlFJ7LjD0EBH1n8MbQnm1A5+dakU4JqMoMzUXwyHqCcMP0SgUjknYfaIF5VWO+J49qTQ2u3NzUu250HP91NyUXNCBiKg/2sMx7KxtRsWxjuHR+RYDLClUO4j6guGHaBQRQuCYw4cPD9lxzO6FWa/BtLzU2bOnM/Rozm1OumgaQw8R0aXEJBkH6t0oq3bgVHMAGQZNyg2PJuorhh+iPpJkgT11rXD6QrCa9bi6OGtI90Ro8YdRVuXAp3UtiMQEJuakpczyo53D27QMPUREfRaTZPzv/kZU1Djh8oVhTddhSm4a1Cm2AihRfzD8EPXBpsNNWLOxCk2eUPxYvkWPp2+fjltK8wf12JGYjD11rSivtqPRHUK+RY8Mo3awTR4T/KEYmjxBaNRKzBn/xfA2JTfaIyK6qLf2nMHPP6iGNxSLHzPpVOfmRpoT2DKi0Y3hh+gSNh1uwg/X74O44LjdE8IP1+/DunvnDCgACSFwwtWODw83obrRC4NWhWk2c0p88D8/9Mxm6CEi6jN/OIbfbanB7z+q6+E2Ce8fsmP5TDAAEfWC4YfoIiRZYM3Gqm7BBwAEAAWANRursHR6Xr+GwHkCUWw56sDHtc0IRiQUZaXGSjznh54rx2fg+qm5mGpNjcBHRDQYUUnGvtNtKKt24I3dpy96bkWNixs/E/WC4YfoIvbUtXYZ6nYhAaDJE8KeulbMn5x9yceLSTI+O92GzUfsaGgLwpauR2GGIekLlD8cQ5M7CI2KoYeIqD86F8Ipq3LgaJMP7mAE4Zh80fv4wxLOuoMYl2kcoVYSjR0MP0QX4fT1Hnz6e96p5o4hbofOeqBTqzDVZh7SBRNGo/NDzxVFGVg0jaGHiKivmjxBlFc5UHm6DRFJxrhMIyT54sGnU3s4dumTiFIQww/RRVjN+kGf5wtFse2oCzuOu+APxzAu0wCjNrn/1zs/9MwqysCiqbkpM5+JiGiwvKEodtS4sON4M9oCkS779aTp+lY/+noeUarp11qI69atw6xZs5Ceno709HTMnz8fH374Yfz2++67DwqFosufL33pS10eIxwO4+GHH0ZOTg7S0tJwxx13oKGhYWieDdEQu7o4C/kWPXr7yK5Ax6pvVxdndbtNkgUqT7fiN1uO471DjdCqlSixmpI6+LSHYzju8MHlC2FWUQb+4YbJ+D/XTcLl+ekMPkRElxCJydhV24xfldXgnYONUACYZjN32ai0MMMAk+7ic0RNOhUKMwzD3Fqisalfn8LGjRuHX/ziF5gyZQoA4A9/+AO+/OUvY//+/ZgxYwYA4JZbbsFrr70Wv49W23XJ3lWrVmHjxo3YsGEDsrOzsXr1aqxYsQKVlZVQqZJ/wjeNLSqlAk/fPh0/XL8PCqDLwgedH+Wfvn16t+FrDW0BbDpsx/4zbqiVCpTkmpJ634VAJIZGdwgqJdjTQ0TUT7IsUNXkRVmVAzUOH/QaFUqsJqiV3euGQqHAoqm5eP+QvdfHWzQ1N+nnkhINVL/Cz+23397l+5///OdYt24dPvnkk3j40el0yMvL6/H+Ho8Hr776Kt544w0sWbIEALB+/XoUFRWhvLwcN99880CeA9GwuqU0H+vundNtn5+8Hvb5aQ/H8FGNC9uPOeEJRTEuw5jUQw/CUQkN7iAAYHp+Om64LBeX57GXh4ior+pbAyivduDAGTdissCELCN0l1j9c4rVjOUzO1Z184el+HHu80N0aQP+VCZJEv785z+jvb0d8+fPjx/fvn07rFYrMjIysGjRIvz85z+H1WoFAFRWViIajWLZsmXx8wsKClBaWopdu3Yx/NCodUtpPpZOz8OeulY4fSFYzR1D3Tp7fGRZ4POzHvztSBNOutqRZdRiqtWctFfeIjEZje4gorJAidWEGy+zorQgPal7t4iIhpI7EMH2Yy7sOtEMTzCKwgwDzHrNpe94zhSrGZNzTTjrDqI9HEOaTp0Sq4cSDVa/w8+hQ4cwf/58hEIhmEwm/PWvf8X06dMBALfeeiu+/vWvY8KECairq8O//uu/4qabbkJlZSV0Oh3sdju0Wi0yMzO7PKbNZoPd3nv3bTgcRjgcjn/v9Xr722yiQVMpFT0uZ93kCeJvh+2oPN0GhQKYnGuCJklDQEyS0egJIRiVUJxjxI3TbLiyKANadXI+X6KLYW2igQhFJXxa14ptR51odAeRa9Zhmm1gF8sUCgWXsybqp36Hn2nTpuHAgQNwu934y1/+gm9/+9uoqKjA9OnTcffdd8fPKy0txbx58zBhwgS8//77uOuuu3p9TCHERf+nX7t2LdasWdPfphINq1BUwo7jLmw76kJLe7jfV+3GEkkWsHtD8IViKMo04IZpVsybmJkSG7MS9Ya1ifpDlgUOnfWgrNqBWocfJp06JbY8IBpt+h1+tFptfMGDefPmYe/evfj1r3+N//f//l+3c/Pz8zFhwgQcP34cAJCXl4dIJIK2trYuvT9OpxMLFizo9Wc++eSTePTRR+Pfe71eFBUV9bfpRENCCIEjjV5sOmJHrcMHi0E74Kt2o50sBJy+MNyBCPLS9bhtZh6uLs6GKYnnMRH1FWsT9dWp5naUVTnweYMbAFCck8Yec6IEGfQnGCFEl27/87W0tKC+vh75+R0TwufOnQuNRoOysjKsXLkSANDU1ITDhw/j+eef7/Vn6HQ66HS6wTaVaNCc3hA2Vzmw91QrJFmgOMeUlAVMCIFmfwTN/jByzTp8+YpCXDslBxZjcvZsEQ0EaxNdSos/jG3HnPjkRAvaIxIKMwxJvQgO0VjQr/8Df/SjH+HWW29FUVERfD4fNmzYgO3bt2PTpk3w+/145pln8NWvfhX5+fk4deoUfvSjHyEnJwd33nknAMBiseD+++/H6tWrkZ2djaysLDz22GOYOXNmfPU3otHIHYhg5/Fm7KxtRmt7BIUZBqQbki8ICCHgDkRh94aQadTgtpn5uHZKDnLN/IBHRNRXwYiE3Sebse2oEw5vGLZ0PQq4GAHRqNCv8ONwOPB3f/d3aGpqgsViwaxZs7Bp0yYsXboUwWAQhw4dwn/913/B7XYjPz8fN954I9566y2YzV8sufjLX/4SarUaK1euRDAYxOLFi/H6669zjx8alQKRGD492Yptx5xweEPIMmpxWV5yDnHzBqNo8gZh0qpx02VWLJqWi3wLN8kjIuorSRY4UN+G8moHTrraka7XcM8zolFGIYQQlz5tdPF6vbBYLPB4PEhPT090cygJRWIy9p1pw9ajDpxpCcCk1yAvXZ+UE1P94Rga3UHoNUpcWZSJG6blYnyWMSkDHg0Nvgf3jK9L6hJC4ITLj7IqBw6f9UKtUqAww5C0K38SDYfjDh/uuWYCFpbk9Pu+/Xn/5cBTovNIssDhsx5sqXag1umHVq3E5FxTUu5fE4xIOOsOQqUEZo/PwI3TrJhiNTH0EBH1g9MbwtajTuypa0UoKqEw0wCjlh+viEYr/t9JhI6rdsedfmypduBIY8deHeOzjdCpk284Zjgm4WxbEALAZflm3HSZFZfnpXNYBhFRP/jDMXxc24yKYy60+MOwWfQoyuKeO0SjHcMPpbz61gC2HnVi/5k2RGIyCpJ0NZ6oJOOsO4hoTMYUqwk3XmbFzEJLUvZqERENl6gkY9/pjnk9Z1oDyDBoMTXPDCV7zYnGhOT7hEfUR05fCB/VuPDpyVZ4Q1EUWJJzBbeYJKPJE0IgImFijhE3TrNi9vjMpFyim4houAghUNXkxZZqJ47afdCqFJiSpMOiiZIZww+lHE8wio/PLVvd4g/Dmq5HvkWfdHNdJFnA4Q3BG4phXKYed84pxLwJWTBok28oHxHRcKpvDaC82oEDZ9yIyQJFmQboNXwvJRqLGH4oZQQjEj6ta8H2Yy40eYLINGoxLQmXrZaFgMsXRmt7BHkWPW6eYcOXJufAlIRD+YiIhlOLP4yKGhd2n2iBL9wxQsCsT74RAkSphJ+GKOlFYjL2n2nD1mNOnG4OwKRTo8RqTrplq4UQaGmPwOULI8ekw5evLMC1U3KQYdQmumlERGNKIBLD7hMt2H6sY5NSq1mHfEvyXSwjSkUMP5S0JFngSKMHW6qdOO7wQatWYlJuWtLtuyCEgDsYhcMTgsWowa2leVhYkotcsy7RTSMiGlOikoz9Z9zYctSBU83tsBi4SSlRsmH4oaQjhECt048tR504ctYDAWB8lhG6JByf7Q1G0eQJIU2nwo2XWXH91FwUZBgS3SwiojFFCIHqJh+2HHXgaJMPapUCk3NNSXexjIgYfijJ1LcGsP2YE/vOtCEcTd5lq0NRCQ1tQahVCnxpUhZumGbFhGwjh2QQEfVT53YH+860ISYJjONiBkRJLfk+FVJKcvnC2HG8Y1KqNxRDvkUPS1byTUqNyTIa3SGEoxIuy0/H0uk2XJaEizYQEQ231vbIucUMmuEJRlGQYUA6FzMgSnoMPzSmeUNR7K5tRkVNM5r9HZNSp9lMSRcGROcKboEIijKNWHy5DXMncK8eIqL+CkRi+PRkK7Ydc8LuCSHXrMM0Gy8iEaUKhh8ak4IRCXtPdRSvRncQGeeWrU7GHba9wSgazy3NfccVBbhuai6vThIR9VNMkrG/3o0t1Q6cammHWcfFDIhSEcMPjSlRScaBeje2HnWgzhVAmk6VlMtWA1/M69GoFFg4JQeLL7dxMQMion4SQuCYw4fyKieqmzxQK5WYlMPFDIhSFcMPjQmyLFDV5MWWagdqHH6olYqkXLYa4LweIqKh0tB2bjGD022ISDLGZRph4GIGRCmN4YdGNSEETrjase2oE583uCEEknYlHs7rISIaGm3tEXx03IWPa5vhCcZQYNEj3cDhwkTE8EOjlBACZ1oDqKhx4cAZN4JRCYVJumw10HVez5evKMTCqTmc10NE1E/BiIRP6lqw7WjHYgY5puRcBIeIBi45P0nSmHbWHcRHx1z47HQr/OEY8i0GFGUZE92sYXHhvJ4l023It3BeDxFRf8QkGQcb3Civ7pgPatKrMdWWnPNBiWhwGH5o1HB4Q9hx3IU9da3wBKOwpetRmGFIyit2nNdDRDR4QgjUOPwor3aguskLVRLPByWiocHwQwnX7A/j49pm7D7Rgrb2SHzPBaCjF6g9HEOaTp0UQYjzeoiIhkajO4gt1Q7sO9OGcPTcYgba4Z8PKoRIutpElEoYfihh2toj2HWiGR/XtqDZH+4Ym32u96PW6UNFjQv+sBQ/36RTYdHUXEyxmhPY6oHjvB4iosFzByL4qKYZH9e64A5GkW8xwJI1Mu+lyVibiFINww+NOE8wik9PtmDH8WY4vCFkpXXdoLTW6cP7h+zd7ucPS3j/kB3LZ2JMFRnO6yEiGrxQVMKnda3YdrRjc+uOxQxGbrhwstUmolTF8EMjxh+OYW9dKypqXGjyBGExdN9dWwiBihrXRR+nosaFybmjf/UezushIhq8QCSGA2fc2FHr6ljMQDfyixkkU20iSnUMPzTsghEJn51uRcUxF+rbAjDrNCix9ly4zrqDXYYT9MQflnDWHcS4zNG5Ahzn9RARDZ4vFEXl6TbsPN6Ms+4gdGolinPSEvJemgy1iYg6MPzQsAlFJew704aKGhdONwdg1KowJdcE9UVW4WkPx/r02H09b6RxXg8R0eC4AxHsqWvFrhMtaPIEkaZVozgnsSu4jfXaRERfYPihIReJdey3sP2YEydd7dCplX1eerSvm5iOts1OOa+HiGhwXL4w9tS1YPeJFrj8YaTrNZhiNUGtTHyv+VitTUTUHf8vpSETk2QcOuvB9mMuHHf6oVYqMDG7f0MUCjMMMOlUFx1eYNKpUJgxOoIF5/UQEQ1OkyeIT0604NO6VrS2R5CZpu11aHSijLXaRES9Y/ihQZNkgapGL7bXOHG0yQelAhifaYBO0//9FhQKBRZNze1xRZ1Oi6bmJjxccF4PEdHACSFQ3xrErhPN2HemDe5ANL7dgXIUXjwaK7WJiC6N4YcGTJYFjjl82H7MiapGL2TRcXVssJvMTbGasXwmRu1eChfO67luag7MnNdDRHRJQgicbG7Hx7XNOFjvhj8cQ65JNyZ6zEd7bSKivmH4oX4TQqDW6cf2Yy4cPutBVJZRYDEM6VjnKVYzJueaRtUu2uGYhPpWzushIuovWRaocfqw83gzjjR6EIzKsJl1CX9f76/RWJuIqH8YfqjPhBA41RJAxTEnDja4EYp2hB6Tfnj+GSkUilGxZKgsCzR5QmiPxHBZnhnLZuSNiauURESJFpNkVDf5sLPWhaNNPkQlGbZ0PcZnjd3e8tFSm4hoYBh+qE/qWwOoqHFh35k2BMIS8i16pI/h4tUXQgi4A1HYvSHkW/T48uwCXF2cBZ16cMP6iIiSXSTWsQDOzuMdC+AIAeRb9FwNjYgSju9CdFFNniB2HG/G3rpWeEMx5KXrMC4FuviDUQkNrQHotSosnW7D4sttyErTJrpZRESjWigq4WC9GzuON6Ou2Q+FQoECy+DnghIRDRWGH+qR0xfCzuPN+PRkC9zBKKxmPfIt+qQPPZIscNYdRDgmY0aBBTeX5mFyblrSP28iosFoD8ew/4wbO2pdONMSgEalRFGmcUCrfhIRDSeGH+rC6Q3h05Ot2HWiGa3tEeSYdZhmS/75LUIItLRH4PKFMS7TgKXT8zBvYmZCdxQnIhrtPMEo9p1uw47jLjS6QzBoVf3e342IaCQx/BAAoNEdxO4Tzdh7qg2t7RFkp2kxLUUm9beHYzjrDsKkU2P5rHzcMM0KiyG55zMREQ1Ga3sEe+ta8PGJFji8IZh0akzOTYOaF4yIaJRj+ElhF24y5wlEkT1G9lsYCjFJRkNbEJIQmD0+AzfPyMOE7LREN4uIaNRyekPYU9eK3Sdb4PKFkWHUoMRqhkqZ/DWDiJIDw08KEkKg7twmcwfObTJnNesxLS/55/QAHc/f5Q+jtT2CCdlGLJueh9njM1m8iYh6cdYdxCcnmrHnVBva2iPIStNims0MJd83iWiMYfhJIbIsUOvyY2dtMw43eBCISLClj71N5gbDH4qhwR1EhlGDL19RiOun5cLEpVeJiLoRQuB0SwC7TjRj/xk3vKEoclJodAARJSd+6ksBkixw1O7FzuPNqGryIhKTkZeux7jM1Ak9UUlGfVsAEMCXJmVh6XQbN6kjIuqBJxDFUbsXB+rdOObwoT0swWrWpcSKn0SU/Bh+klhMklHV5MWOmmYctXshyQL5FgNM+tT5tQsh4PCG4QlGMNlqwrLpeZhZaOFQDSKi84SiEmocPhw668Hhsx60tUehVimQY0qNvd2IKHWkzqfgFNK5s/aO4y7UpvDO2t5gFI3uILJNOnx17jgsnJLLjfaIiM6JSTJOtbTjSKMX+8+44fSFIYRAplGLKVYT50ESUVLq15qU69atw6xZs5Ceno709HTMnz8fH374Yfx2IQSeeeYZFBQUwGAw4IYbbsCRI0e6PEY4HMbDDz+MnJwcpKWl4Y477kBDQ8PQPJsUF4pK+PRkC36z9Tj+c+dJnHD6UZBhwBSrKaWCTzgmodbpR2sggoUlOXhkSQmWTs9j8CGilCdEx0bOW6odeLGsBr/dUosPDjUhEIlhQpYRU21m5Jp1DD5ElLT69Yl43Lhx+MUvfoEpU6YAAP7whz/gy1/+Mvbv348ZM2bg+eefx4svvojXX38dU6dOxc9+9jMsXboUx44dg9lsBgCsWrUKGzduxIYNG5CdnY3Vq1djxYoVqKyshErFD6cDEYxI2H+mDR8d79hZW52iO2vLskCTN4T2cAwlNhNunpGH6fnpHK5BRCmvtT2Co00dPTwnm/3wh2Iw6tSwputg1KbOxTEiIoUQQgzmAbKysvDv//7v+O53v4uCggKsWrUKTzzxBICOXh6bzYbnnnsOP/jBD+DxeJCbm4s33ngDd999NwCgsbERRUVF+OCDD3DzzTf36Wd6vV5YLBZ4PB6kp6cPpvljmj8cQ+XpNuyocaGhLQi9Rol8iyEld9ZuC0Tg8IRgs+ix+HIbrinOgj7Fwh/RSOF7cM9G2+sSiMRwzN4xj6eq0Qt3IAKNSokckw5mvZoXhohoVDnu8OGeayZgYUlOv+/bn/ffAV/ukSQJf/7zn9He3o758+ejrq4Odrsdy5Yti5+j0+mwaNEi7Nq1Cz/4wQ9QWVmJaDTa5ZyCggKUlpZi165dvYafcDiMcDjc5QmmMk8wis9OtWLH8WbYPUEYtWpMyk2DJgV31g5FJdS3BqDXqHDjZVYsmW5DjkmX6GYRUQoYjbUpKsk46WrHkUYPDtS74fKHoRBAVpoWJVbuy0NE1O/wc+jQIcyfPx+hUAgmkwl//etfMX36dOzatQsAYLPZupxvs9lw+vRpAIDdbodWq0VmZma3c+x2e68/c+3atVizZk1/m5p0Wtsj2FvXgo9PtMDhDcGs12Cy1QS1MvVCjyQLNLqDCMUkXF6Qjptn5KHEauKVTCIaMaOlNslyxzyeqkYP9p1xo9EdRFSWYTFoUZydmhfGiIh60+/wM23aNBw4cAButxt/+ctf8O1vfxsVFRXx2y/88CmEuOQH0kud8+STT+LRRx+Nf+/1elFUVNTfpo9ZTl8Ie062YvfJFjT7w7AYNCixmlNyQqoQAq3tETh9YRRmGPC16eMwb2JWSg71I6LESnRtcvnCOGrvmMdzqrkd7ZEYTDo1CjIMHPZLRNSLfocfrVYbX/Bg3rx52Lt3L37961/H5/nY7Xbk5+fHz3c6nfHeoLy8PEQiEbS1tXXp/XE6nViwYEGvP1On00GnS72hTE2eIHafaMGeula0tkeQlabF1BQethCIxNDQFkSaVoVbSvNw02VWZBi1iW4WEaWoRNQmXyiKGocPB850bEDqDcagVSuRY9Km1MbVREQDNeglXoQQCIfDKC4uRl5eHsrKyjB79mwAQCQSQUVFBZ577jkAwNy5c6HRaFBWVoaVK1cCAJqamnD48GE8//zzg21KUhBCoKEtiF0nWlB5uhWeQBTZJh0uyzOnbFGLyTLOtgURlQRmjbPg5hl5mJRrSnSziIhGROfy/Ucavfi8wY0WfwRKhQLZaVrYbHooU7Q2EBENRL/Cz49+9CPceuutKCoqgs/nw4YNG7B9+3Zs2rQJCoUCq1atwrPPPouSkhKUlJTg2WefhdFoxD333AMAsFgsuP/++7F69WpkZ2cjKysLjz32GGbOnIklS5YMyxMcK4QQONUSwMe1zThwpg2+cAy5Jh2mpXDoiQ9x84cxIcuIpdPzMGd8BtQcv05ESU6WBU63BlDd5MG+0240eUKQhUCGQYNJOWl8HyQiGqB+hR+Hw4G/+7u/Q1NTEywWC2bNmoVNmzZh6dKlAIDHH38cwWAQDzzwANra2nDNNddg8+bN8T1+AOCXv/wl1Go1Vq5ciWAwiMWLF+P1119P2T1+ZFmg1uXHrtpmfN7gQSAag9WsR0FGag9fCMcknGkJwKBV4bbSfNx0uRXpek2im0VENKycvhCqGr2oPN2G+tYAghEJZoMGRVkG6NSpWSeJiIbSoPf5SYTRtpfCQEiyQHWTFztrm1Hd5EUkJiMvXZ/yey8IIeDwheEJRHFZvhkrZuVjitV86TsSnSPJAnvqWuH0hWA163F1cVZKLg4ynJLhPXg4DPZ1CUYk/N/Nx9DQFoBBo0KOSQejVpXSNYEoWQjRsSpjeziGNJ0ahSl+kbsno36fHxqYSEzGobMefFzbjOMOH2QB5KXrYdLzVxGIxFDfGkSGUYOvzi3EdSW5XLGI+mXT4Sas2ViFJk8ofizfosfTt0/HLaX5F7knUeLFZBnt4Rjy0vVczIUoidQ6faioccEfluLHTDoVFk3N5QXeBOAn7hESikrYf8aNnbUunGpuh0KhQL5FD6OWvwJZdOzZE4xIuHJ8BpbPzEdRljHRzaIxZtPhJvxw/T5c2JVt94Tww/X7sO7eOQxAREQ0omqdPrx/qPtelv6whPcP2bF8JhiARhg/eQ8zXyiKfWfc2HnchYa2ILQqJYqyjBy7fY4vFMVZdxC2dD3umjMO1xRncSIv9ZskC6zZWNUt+ACAAKAAsGZjFZZOz+MQOCIiGhFCCFTUuC56TkWNC5NzuUn7SGL4GSZt7RHsPdWKXSdaYPcEYdSqUZzDnbY7xWQZDa1ByELg2ik5uLU0H7nm1NvLiYbGnrrWLkPdLiQANHlC2FPXivmTs0euYURElLLOuoNdhrr1xB+WcNYdxLhMjngZKQw/Q8zp7fiAtftkC5r9YaTrNZhsNUGtZOjp1BaIwO4JYUK2EbfOzMeV4zJSduNWGhpOX+/BZyDnERERDVZ7ODak59HQYPgZAp0reHxysgV7T7WhrT2CrDQtplrN/FB/nqgk40xLABq1Estm2LB0eh4sBi5fTYNnNeuH9DwiIqLBStP17WN2X8+jocFXexCEEDjZ3I7dJ1pwoN4NXzCKHLMOl6XwxqQ9EUKg2R9BS3sEJTYTls/M52tEQ+rq4izkW/Swe0I9zvtRAMizdCx7TURENBIKMwww6VQXHfpm0qlQmGEYwVYRw88AyLJAjdOHncebcaTRi2BEgjVdhwKLnh/oLxCKSjjTGoBZr8aXryjADZflcoU7GnIqpQJP3z4dP1y/DwqgSwDq/D/y6dunc7EDIiIaMQqFAoum5va42lunRVNz+dlxhPFTaD/EJBlVTV7sPN6Mo3YvYpKALV2P8VyWuRshBJo8IfjDMZQWWLB8Vj4m5qQlulmUxG4pzce6e+d02+cnj/v8EBFRgkyxmrF8JrjPzyjC8NMH4ZiEzxs82Hm8GSdcfgjRsXEix2j2rD0cQ31bEDkmLe6eVYQFk3OgVXPBBxp+t5TmY+n0POypa4XTF4LV3DHUjT0+RESUKFOsZkzONeGsO4j2cAxpOjUKMwzs8UkQfnq/iEAkhgNn3NhR68Lp5gDUKgUKMwzQa7hHT08kuWPhh6gk45riTNw6Mx/5Fo5jpZGlUiq4nDUREY0qCoWCy1mPEgw/PfAEo9h3pg0fH29GQ1sQeo0S47O5MenFeIJRNLqDKMw04NbSfMydkMmr7UREREQ0qjD8nKfFH45vTOrwhmHSqTA5Nw1qbkzaq5gk40xbAEoocOM0K24uzUNWmjbRzSIiIiIi6obhB4DdE8InJ1vwaV0LWtsjsBg0KLGa2HNxEUIItLZH4PKHUZyThuUzC1BamM7xq0REREQ0aqVs+BFCoL41iE/qWvDZqVZ4AlFkmbSYajNDyQ/wFxWOSTjTEoBBq8JtM/Ox+HIbTFz8gYiIiIhGuZT9xPr+503YdswJfziGXJMO07jp5iUJIeDwheEJRHFZvhkrZuVziUYiIiIiGjNSNvw0uIMIx2RMszH09EUgEkN9axAZRg2+OrcQ15XkctU7IiIiIhpTUjb8AIBaqWDwuQRZCDS6gwhGJFw5PgPLZ+ajiJu6EhEREdEYlNLhhy7OF4rirDsIW7oed80Zh2uKs7jyHRERERGNWQw/1E1MltHQGoQsBK6dkoNbS/ORa9YlullERERERIPC8ENddG5WOiHbiFtn5uPKcRlQcslvIiIiIkoCDD8EAJBkgfq2ACCAmy6z4tbSfFiMmkQ3i4iIiIhoyDD8UJfenhWzCjBrnIULQRARERFR0mH4SWGdvT1CCPb2EBEREVHSY/hJUd5gx0pu47OMWD4rH1cWZbC3h4iIiIiSGsNPipFkgfrWAGQhsGhaLpbPzEeGUZvoZhERERERDTuGnxRyfm/PbbPyMZu9PURERESUQhh+UkDn3B5ZANdP7ejtyUxjbw8RERERpRaGnyTnC0Vxti2IwkwDls8qwOwi7ttDRERERKmJ4SdJSbJAQ1sAkiywsCQHy2cVIIu9PURERESUwhh+kpAv1DG3pzDDgNtm5mPO+Ez29hARERFRymP4STAhBM66g2gPx5CmU6MwwzDgRQjkc3N7YrLAtVNysHxmPrJNuiFuMRERJbuhrE1ERKMJw08C1Tp9qKhxwR+W4sdMOhUWTc3FFKu5X4/lD8XQ0BZAQYYBt83Kx1z29hAR0QAMZW0iIhptlIluQKqqdfrw/iF7l+ICAP6whPcP2VHr9PXpcWRZ4ExrAA5fCAum5ODBm6bgqolZDD5ERNRvQ1WbiIhGK/b8JIAQAhU1roueU1HjwuRc00WHGfhDMdR39vbMzMe8CeztISKigRmq2kRENJox/CTAWXew21W1C/nDEs66gxiXaex2mywLNLiDiEoyFkzJxvKZBcg1c24PEREN3GBrExHRWMDwkwDt4diAz2sPx3CmNYCCDD1uLc3HvIlZULG3h4iIBmkwtYmIaKxg+EmANF3fXvbzz5OFwNm2ICKSjPmTs7FiFnt7iIho6AykNhERjTV8B0uAwgwDTDrVRYcXmHQqFGYYAHRcZatvDcBm6ejtubqYvT1ERDS0+lubiIjGIq72lgAKhQKLpuZe9JxFU3MhADS0BdDkCeGaSVl46KYpmD85m8GHiIiGXF9rExc7IKKxjD0/CTLFasbymeh1L4V8iwE1Dh9s6Xp8dW4erikemdAjyQJ76lrh9IVgNevZy0RElEIuVZsSsc8PN1wloqHUr/Czdu1avP322zh69CgMBgMWLFiA5557DtOmTYufc9999+EPf/hDl/tdc801+OSTT+Lfh8NhPPbYY/jTn/6EYDCIxYsX46WXXsK4ceMG+XQurfPD/ZFGD1rbIxifZUzYm+gUqxmTc01d3tTzLXo0eUJo8oRwdXEWVswqgC1dPyLt2XS4CWs2VqHJE4ofy7fo8fTt03FLaf6ItIGIKBV11qZTLe2w6DWwGDSjqjYlKnBww1UiGmr9Cj8VFRV48MEHcdVVVyEWi+Gpp57CsmXLUFVVhbS0tPh5t9xyC1577bX491qttsvjrFq1Chs3bsSGDRuQnZ2N1atXY8WKFaisrIRKpRrkU+pdTx/uD9a7E/omqlAo4kuGtodjOO70w2rW46tz8nB1cRbUqpEZmbjpcBN+uH4fxAXH7Z4Qfrh+H9bdO4cBiIhoGPRUmz4+0TxqalOidG64eqHODVeXzwQDEBH1W7/Cz6ZNm7p8/9prr8FqtaKyshLXX399/LhOp0NeXl6Pj+HxePDqq6/ijTfewJIlSwAA69evR1FREcrLy3HzzTf39zn0re29fLgfDW+iQgg0ukMIRCXMm5CF268oQJ5lZHp7gI4rjms2VnV7bQBAAFAAWLOxCkun53EIHBHREBrNtSmRuOEqEQ2XQXUreDweAEBWVlaX49u3b4fVasXUqVPx/e9/H06nM35bZWUlotEoli1bFj9WUFCA0tJS7Nq1q8efEw6H4fV6u/zpj4t9uO9UUeOCEBc7Y3gEIxKOOfzQa5S45+rx+M61E0c0+ADAnrrWLlccLyQANHlC2FPXOnKNIiIa5ZK5NiVafzZcJSLqjwGHHyEEHn30USxcuBClpaXx47feeiv++Mc/YuvWrXjhhRewd+9e3HTTTQiHwwAAu90OrVaLzMzMLo9ns9lgt3fv3gY65hpZLJb4n6Kion619VIf7oGRfxMVQqDJE0RDWwBzxmfgwZumYGFJzogNczuf03fx16a/5xERpYJkrE2jBTdcJaLhMuBP2g899BA+//xz/OlPf+py/O6778by5ctRWlqK22+/HR9++CFqamrw/vvvX/TxhBC9dl0/+eST8Hg88T/19fX9amtfP7SP1JtoKCrhuNMPlUKBlVcV4bsLi5FvSdy+CVZz33qa+noeEVEqSLbaNJpww1UiGi4Detd4+OGH8e677+Kjjz665Apt+fn5mDBhAo4fPw4AyMvLQyQSQVtbW5feH6fTiQULFvT4GDqdDjqdbiBNBdD3D+3D/SYqhIDTF4Y7EMXMQgvuuLIARVmJnVAKAFcXZyHfoofdE+px+IUCQJ6lY9lrIiLqkCy1aTTihqtENFz61fMjhMBDDz2Et99+G1u3bkVxcfEl79PS0oL6+nrk53esFDZ37lxoNBqUlZXFz2lqasLhw4d7DT+D1fnh/mJTIof7TTQc6+jtkYXAXXMK8b3ri0dF8AEAlVKBp2+fDgDdXqPO75++fToXOyAiGkKjoTaNVtxwlYiGS7/Cz4MPPoj169fjzTffhNlsht1uh91uRzDYMR7Z7/fjsccew+7du3Hq1Cls374dt99+O3JycnDnnXcCACwWC+6//36sXr0aW7Zswf79+3Hvvfdi5syZ8dXfhtrFPtx3Gq43USEEXL4w6prbcVmeGT+8YQqWzciDTj18S3oPxC2l+Vh375xuiy3kWfRc5pqIaBgksjaNBR0brubBpOtaL006FZbPzEvJVfCIaPAUoh/LyPT2Bvzaa6/hvvvuQzAYxFe+8hXs378fbrcb+fn5uPHGG/Fv//ZvXSaChkIh/PM//zPefPPNLpuc9nWyqNfrhcVigcfjQXp6el+b3+NeCsO5WVpUknG6JQCDRokbL7Ni8eU26DWjK/RcqHOjPacvBKu5Y6gbe3yI6HwDfQ9OdmOlNo01QohRseEqEQ2v4w4f7rlmAhaW5PT7vv15/+1X+BktBlN4Oz/cv7nnNFrbI5g7PnNY3kRb2yNw+kIosZlxxxUFmGpjASOi5MDw07PB1qbtx5x4ZcdJWPQaTMsz8wM+EaWUkQo/KTeLUqVUYP7kbBxscONok3fIi0tMknGmNQC1SolbSvOxbLotJSerEhFR36mUClxdnIWyKgd0aiWDDxHRMOGn8iHkDkTQ5AlhUm4abr+iANPz01nAiIiIiIhGCYafISDJAvVtAQDAkuk23FKah3S9JsGtIiIiIiKi8zH8DJI3GMVZdxATso1YMasAs8ZZ2NtDRERERDQKMfwMkHyut0eSBRZNy8XymfnIMGoT3SwiIiIiIuoFw88A+MMxNLQGUJBhwPJZ+ZgzPhNKLgdNRERERDSqMfz0gywEzrYFEZFkzJ+SjRWzCpBj0iW6WURERERE1AcMP30UiMRwpjUAW7oet5bmc/NPIiIiIqIxhuHnEoQQaPKE0B6OYd6ELNxxZQFs6fpEN4uIiIiIiPqJ4eciQlEJp1rakZ2mwzeuHo8Fk7OhVikT3SwiIiIiIhoAhp8eCCHg8IbhCUVxRVEGvnxlIQozDIluFhERERERDQLDzwXCMQmnWwJIN2jw9bnjcF1JLrRq9vYQEREREY11DD/nCCHQ7I+gtT2M6QUWfPnKAkzITkt0s4iIiIiIaIgw/ACIxGScbmlHmk6Nr8wuxA3TrNBrVIluFhERERERDaGUDz8t/jBc/jCm5ZlxxxWFmGI1JbpJREREREQ0DFI6/DT7I1AqFVg+swBLplth1Kb0y0FERERElNRS9tO+QaPCrHEW3HFlAabZzFAouGEpEREREVEyS9nw85XZhdCoFOztISIiIiJKESn7yd9i0CS6CURERERENIK4gQ0REREREaUEhh8iIiIiIkoJKTvsjYiIiIiI+k4IAQFAiI6/y51f0fV7IQD5vK+yAATEuftdeFvHY0YkeUSeA8MPEREREdEYIi4IFrIQkGTR5XtZ7vi7FP/7udvO/b3jOACIS656LDr/owAUABQKQKlQQIFzXxXdvyoUCigVgALnviqUUCkVUCsVUKkU0CiVUKs6vteolNCqlCjMNAz3S8fwQ0REREQ0VM4PJlI8lHSEDum8UNLZ6yHL5wLKBbeJ88JG18fv+NoZNFTKjq9KhQJKZUfQUJ37u0oJ6FUqaNQd4UKnVkKnVkGr7vi7Vq2EWqWEWqk491gdYUSpVEB17nuV8oufoVYqoVTii/O6nKOI/9weH+PcbYnG8ENEREREKaczbEjnB5ALAkpnb4p0QW9JZzi5sMPkXC6BAh0BQano/KqIf+0MLFplRwjRdgkmHd/rNR23qZVKaFQd99WcCylqVcexzts6w4tG1dGTEu9ROXeOahQEjtGE4YeIiIiIRj1xLoR09qZ0+9PDcbm33hOIc8OxvujV6OyhOD+g6NQq6DQdvSU6tRI6TUcw0Z87rlYqoVV3BpOOv3eEkvP+rlZCcy6caNRfBBXluaFhNLIYfoiIiIhoWHROgu81tPQSZHqahyLQMZxL3RlSzhtOpVYpYdIoYdB2BBODVgWjVgWDVg2t6lzwONczolUp470nvf2983z2miQfhh8iIiIi6lXnsK+YJBCTZcTkrt9Lsogfu5AA4nNQuswNOfdHp1RBr1XBoDkXWM591WtU0KqUXwwLOzc0rMvX84+rlKNiPgmNfgw/RERERElOiI6AEpMFYtIXgSUmC0jnhZqYJLrNYwEQ723pnF+iViqQplPBqNXBpFMhTac5970qPkRMc2FQUX0xp6Xze5VSwaFfNKIYfoiIiIjGAPm8Hpbze1zODzCdxy4kgHho6Zw0r1IqYNSokGZSw6hVwaRXw6RVx3teOua4dHzVa1TQn5v70vlVo2JwobGH4YeIiIhoBHRO2P8iwHTthfkizMiQ5e4riQHnAoxKAZVSGe+NMWhUSDOpYNSqYdKpYNJpzgWYrmGl8+v5k/YZYCjVMPwQERER9ZEQAlFJIHoutJwfWM4PMRfOf+n8Lj587LyvGrUS6RoVjDo10rQqpJ0bRta59PGFwaXzWOfftSolAwxRHzH8EBERUcqTZIGIJCMakxGVZESlc99Lcrd5MOfvqdIZYIxaDYy6zt6XjhBjOG/+y/nLJceHlJ0LMJysTzRyGH6IiIgoKXX20nSGmI5g88X357aAAdCxIlnn8sYalRImvQrpBgMyjRpkGDRI06njwcZ43sR+/bnwolYpE/lUiaiPGH6IiIhozOicN9Ml1EgyIueCjSyLeKIRQJc9XrRqJbJMamQYtMgwamDWd65Qpu4YaqZVnws5HT0zHEpGlHwYfoiIiGjU6FySORSVEIrKCMc6vkpChgIKCAGolIiHGY1KCYtBA4tBg0yjFhajBmnajgDTGWQ6emvUMGhU3LSSKMUx/BAREdGIEqKj1yYUlRGKSgjHOkKOLANQdMyp0Z+bG1OQYYAtXY8ckw5mfddQ09lTo1VzyBkR9Q3DDxEREQ05WRYIx2SEYhLC50JORJIBAEIAWrXi3N4xHQEnz6JHlrFjOFqGUYP0c705OrUqwc+EiJIJww8RERENiBQfntbRexOKSohKcnyuTOfKZia9CsW5abCZdchM0yLDqIXF0LGQQLpBw6FoRDRiGH6IiIioR73PvxFQoGMTTv25Fc+y0rSwpethTdfFFxSwGDSwGDUw69RcPICIRgWGHyIiIoIsBIIRCe2RGNrDXwxR623+TWe4yTBqkGHQQq/h6mhENPr1K/ysXbsWb7/9No4ePQqDwYAFCxbgueeew7Rp0+LnCCGwZs0avPzyy2hra8M111yD//iP/8CMGTPi54TDYTz22GP405/+hGAwiMWLF+Oll17CuHHjhu6ZERERUY9kWSAQDzoxRM9t4mnQdKyMNqMgDROyjfGQk27oCDhcWICIxrp+hZ+Kigo8+OCDuOqqqxCLxfDUU09h2bJlqKqqQlpaGgDg+eefx4svvojXX38dU6dOxc9+9jMsXboUx44dg9lsBgCsWrUKGzduxIYNG5CdnY3Vq1djxYoVqKyshErFiY1ERERDRZIFAud6c9ojMcTOBR2jVgWTXo0p1gwUZRlhNevivToMOUSUrBRCCDHQO7tcLlitVlRUVOD666+HEAIFBQVYtWoVnnjiCQAdvTw2mw3PPfccfvCDH8Dj8SA3NxdvvPEG7r77bgBAY2MjioqK8MEHH+Dmm2++5M/1er2wWCzweDxIT08faPOJiGgA+B7cs8G+Lr5QFD9/vxo6tRIZRu2A2hCT5PN6dKT43Jw0rQpmgwbjs4woyjTCmq6D1axHjkkLtYpBh4jGtv68/w5qzo/H4wEAZGVlAQDq6upgt9uxbNmy+Dk6nQ6LFi3Crl278IMf/ACVlZWIRqNdzikoKEBpaSl27drVp/BDRESU6qKdQSccQ3skBlkWUCoVSNOqYTFoMLPQgsJMA6zmjkUIstN0XFWNiFLegMOPEAKPPvooFi5ciNLSUgCA3W4HANhsti7n2mw2nD59On6OVqtFZmZmt3M673+hcDiMcDgc/97r9Q602URERENiJGtTJCajPRJDICwhEJEgQ0ClANJ0amQYNZg9PgMFGV8EnSyjFkoGHSKibgYcfh566CF8/vnn2LlzZ7fbLlztRQhxyRVgLnbO2rVrsWbNmoE2lYiIaMgNV22KSDJa2yNoD8cQjEoQAlCrFEjTqZFj1qI4Jw35lo5V16zmjgUJuMoaEVHfDCj8PPzww3j33Xfx0UcfdVmhLS8vD0BH705+fn78uNPpjPcG5eXlIRKJoK2trUvvj9PpxIIFC3r8eU8++SQeffTR+PderxdFRUUDaToREdGQGOrapFAooFUr0dYegUmnRkGGARNzjMi3GGA1d8zRSTdwvxwiosHoV/gRQuDhhx/GX//6V2zfvh3FxcVdbi8uLkZeXh7Kysowe/ZsAEAkEkFFRQWee+45AMDcuXOh0WhQVlaGlStXAgCamppw+PBhPP/88z3+XJ1OB51O1+8nR0RENFyGujalaVX4xlXjodMoYTXrYOLGoEREQ65f4efBBx/Em2++iXfeeQdmszk+R8discBgMEChUGDVqlV49tlnUVJSgpKSEjz77LMwGo2455574ufef//9WL16NbKzs5GVlYXHHnsMM2fOxJIlS4b+GRIREY0BCoUC0wu4eh4R0XDqV/hZt24dAOCGG27ocvy1117DfffdBwB4/PHHEQwG8cADD8Q3Od28eXN8jx8A+OUvfwm1Wo2VK1fGNzl9/fXXuccPERERERENm0Ht85Mo3GOCiChx+B7cM74uRESJ0Z/3X+5sRkREREREKYHhh4iIiIiIUgLDDxERERERpQSGHyIiIiIiSgkMP0RERERElBIYfoiIiIiIKCUw/BARERERUUpg+CEiIiIiopSgTnQDBqJzX1av15vglhARpZ7O994xuEf2sGJtIiJKjP7UpTEZfnw+HwCgqKgowS0hIkpdPp8PFosl0c0YNVibiIgSqy91SSHG4KU7WZbR2NgIs9kMhULR7/t7vV4UFRWhvr4e6enpw9DCsY2vT+/42vSOr03vku21EULA5/OhoKAASiVHT3dibRo+fG16x9fm4vj69C6ZXpv+1KUx2fOjVCoxbty4QT9Oenr6mP9lDye+Pr3ja9M7vja9S6bXhj0+3bE2DT++Nr3ja3NxfH16lyyvTV/rEi/ZERERERFRSmD4ISIiIiKilJCS4Uen0+Hpp5+GTqdLdFNGJb4+veNr0zu+Nr3ja0N9wX8nveNr0zu+NhfH16d3qfrajMkFD4iIiIiIiPorJXt+iIiIiIgo9TD8EBERERFRSmD4ISIiIiKilMDwQ0REREREKSElw89LL72E4uJi6PV6zJ07Fzt27Eh0kxJu7dq1uOqqq2A2m2G1WvGVr3wFx44dS3SzRqW1a9dCoVBg1apViW7KqHD27Fnce++9yM7OhtFoxJVXXonKyspEN2tUiMVi+PGPf4zi4mIYDAZMmjQJP/3pTyHLcqKbRqMQa1N3rE19x9rUFWtTz1iXUjD8vPXWW1i1ahWeeuop7N+/H9dddx1uvfVWnDlzJtFNS6iKigo8+OCD+OSTT1BWVoZYLIZly5ahvb090U0bVfbu3YuXX34Zs2bNSnRTRoW2tjZce+210Gg0+PDDD1FVVYUXXngBGRkZiW7aqPDcc8/h97//PX73u9+huroazz//PP793/8dv/3tbxPdNBplWJt6xtrUN6xNXbE29Y51KQWXur7mmmswZ84crFu3Ln7s8ssvx1e+8hWsXbs2gS0bXVwuF6xWKyoqKnD99dcnujmjgt/vx5w5c/DSSy/hZz/7Ga688kr86le/SnSzEupf/uVf8PHHH/MKdS9WrFgBm82GV199NX7sq1/9KoxGI954440EtoxGG9amvmFt6o61qTvWpt6xLqVYz08kEkFlZSWWLVvW5fiyZcuwa9euBLVqdPJ4PACArKysBLdk9HjwwQexfPlyLFmyJNFNGTXeffddzJs3D1//+tdhtVoxe/ZsvPLKK4lu1qixcOFCbNmyBTU1NQCAgwcPYufOnbjtttsS3DIaTVib+o61qTvWpu5Ym3rHugSoE92AkdTc3AxJkmCz2boct9lssNvtCWrV6COEwKOPPoqFCxeitLQ00c0ZFTZs2IB9+/Zh7969iW7KqHLy5EmsW7cOjz76KH70ox9hz549+Md//EfodDr8/d//faKbl3BPPPEEPB4PLrvsMqhUKkiShJ///Of45je/meim0SjC2tQ3rE3dsTb1jLWpd6xLKRZ+OikUii7fCyG6HUtlDz30ED7//HPs3Lkz0U0ZFerr6/HII49g8+bN0Ov1iW7OqCLLMubNm4dnn30WADB79mwcOXIE69atS/kCA3TM41i/fj3efPNNzJgxAwcOHMCqVatQUFCAb3/724luHo0yrE0Xx9rUFWtT71ibese6lGLhJycnByqVqtuVNKfT2e2KW6p6+OGH8e677+Kjjz7CuHHjEt2cUaGyshJOpxNz586NH5MkCR999BF+97vfIRwOQ6VSJbCFiZOfn4/p06d3OXb55ZfjL3/5S4JaNLr88z//M/7lX/4F3/jGNwAAM2fOxOnTp7F27dqUKTJ0aaxNl8ba1B1rU+9Ym3rHupRic360Wi3mzp2LsrKyLsfLysqwYMGCBLVqdBBC4KGHHsLbb7+NrVu3ori4ONFNGjUWL16MQ4cO4cCBA/E/8+bNw7e+9S0cOHAgZYsLAFx77bXdlp2tqanBhAkTEtSi0SUQCECp7Po2q1KpUmpJUbo01qbesTb1jrWpd6xNvWNdSrGeHwB49NFH8Xd/93eYN28e5s+fj5dffhlnzpzBP/zDPyS6aQn14IMP4s0338Q777wDs9kcvwJpsVhgMBgS3LrEMpvN3caXp6WlITs7O+XHnf/TP/0TFixYgGeffRYrV67Enj178PLLL+Pll19OdNNGhdtvvx0///nPMX78eMyYMQP79+/Hiy++iO9+97uJbhqNMqxNPWNt6h1rU+9Ym3rHugRApKD/+I//EBMmTBBarVbMmTNHVFRUJLpJCQegxz+vvfZaops2Ki1atEg88sgjiW7GqLBx40ZRWloqdDqduOyyy8TLL7+c6CaNGl6vVzzyyCNi/PjxQq/Xi0mTJomnnnpKhMPhRDeNRiHWpu5Ym/qHtekLrE09Y10SIuX2+SEiIiIiotSUUnN+iIiIiIgodTH8EBERERFRSmD4ISIiIiKilMDwQ0REREREKYHhh4iIiIiIUgLDDxERERERpQSGHyIiIiIiSgkMP0RERERElBIYfoiIiIiIKCUw/BARERERUUpg+CEiIiIiopTA8ENERERERCnh/wcTClGd6v+0gAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(10, 4), sharey=True)\n",
    "\n",
    "axes[0].scatter(df_new.iloc[10:][\"Days\"], df_new.iloc[10:][\"Reaction\"])\n",
    "axes[1].scatter(df_new.iloc[:10][\"Days\"], df_new.iloc[:10][\"Reaction\"])\n",
    "\n",
    "axes[0].fill_between(np.arange(10), bounds[0, 10:], bounds[1, 10:], alpha=0.5, color=\"C0\")\n",
    "axes[1].fill_between(np.arange(10), bounds[0, :10], bounds[1, :10], alpha=0.5, color=\"C0\")\n",
    "\n",
    "axes[0].set_title(\"Original participant\")\n",
    "axes[1].set_title(\"New participant\");"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Unnamed: 0</th>\n",
       "      <th>subj</th>\n",
       "      <th>item</th>\n",
       "      <th>site</th>\n",
       "      <th>Y</th>\n",
       "      <th>continuous</th>\n",
       "      <th>dummy</th>\n",
       "      <th>threecats</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.276766</td>\n",
       "      <td>0.929616</td>\n",
       "      <td>0</td>\n",
       "      <td>a</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0.058104</td>\n",
       "      <td>0.008388</td>\n",
       "      <td>0</td>\n",
       "      <td>a</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>-6.847861</td>\n",
       "      <td>0.439645</td>\n",
       "      <td>0</td>\n",
       "      <td>a</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>12.474619</td>\n",
       "      <td>0.596366</td>\n",
       "      <td>0</td>\n",
       "      <td>a</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>-0.426047</td>\n",
       "      <td>0.709510</td>\n",
       "      <td>0</td>\n",
       "      <td>a</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Unnamed: 0 subj  item  site          Y  continuous  dummy threecats\n",
       "0           0    0     0     0   0.276766    0.929616      0         a\n",
       "1           1    1     0     0  -0.058104    0.008388      0         a\n",
       "2           2    2     0     1  -6.847861    0.439645      0         a\n",
       "3           3    3     0     1  12.474619    0.596366      0         a\n",
       "4           4    4     0     2  -0.426047    0.709510      0         a"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = pd.read_csv(\"../../tests/data/crossed_random.csv\")\n",
    "data[\"subj\"] = data[\"subj\"].astype(str)\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "       Formula: Y ~ 0 + threecats + (0 + threecats | subj)\n",
       "        Family: gaussian\n",
       "          Link: mu = identity\n",
       "  Observations: 120\n",
       "        Priors: \n",
       "    target = mu\n",
       "        Common-level effects\n",
       "            threecats ~ Normal(mu: [0. 0. 0.], sigma: [31.1617 31.1617 31.1617])\n",
       "        \n",
       "        Group-level effects\n",
       "            threecats|subj ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: [31.1617 31.1617 31.1617]))\n",
       "        \n",
       "        Auxiliary parameters\n",
       "            sigma ~ HalfStudentT(nu: 4.0, sigma: 5.8759)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula = \"Y ~ 0 + threecats + (0 + threecats | subj)\"\n",
    "model = bmb.Model(formula, data)\n",
    "model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [sigma, threecats, threecats|subj_sigma, threecats|subj_offset]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8a60022ea97640dc99517601576a288e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 4 seconds.\n",
      "There was 1 divergence after tuning. Increase `target_accept` or reparameterize.\n"
     ]
    }
   ],
   "source": [
    "idata = model.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>threecats</th>\n",
       "      <th>subj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>a</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>a</td>\n",
       "      <td>11</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  threecats subj\n",
       "0         a    0\n",
       "1         a   11"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "new_data = pd.DataFrame(\n",
    "    {\n",
    "        \"threecats\": [\"a\", \"a\"],\n",
    "        \"subj\": [\"0\", \"11\"]\n",
    "    }\n",
    ")\n",
    "new_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "p1 = model.predict(idata, data=new_data, inplace=False, sample_new_groups=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlcAAALgCAYAAAC59FJJAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAPTdJREFUeJzt3X9sU/eh//+Xb364SZa4+QE2Vg1kd9nWNqGjoTcjXUe2/OAiKOuld6Gjq9gdd4IB2VyIKBn3qmnV60CmBjbljit6q4WCWKaPtGy9grak2paW5aKlWdEgnbpWTUcy4mXtIjuhuTYL5/tHvz33OgFaJ+/ECTwf0pHqc9523sdy1Wffto8dlmVZAgAAgBF/k+gJAAAAXE+IKwAAAIOIKwAAAIOIKwAAAIOIKwAAAIOIKwAAAIOIKwAAAIOSEz2Bybh8+bIuXLigzMxMORyORE8HAABc5yzL0vDwsLxer/7mb669NjUn4+rChQvy+XyJngYAALjB9PX16ZZbbrnmmDkZV5mZmZLeP8GsrKwEzwYAAFzvwuGwfD6f3SDXMifj6oO3ArOysogrAAAwYz7Kx5H4QDsAAIBBxBUAAIBBxBUAAIBBxBUAAIBBxBUAAIBBxBUAAIBBxBUAAIBBxBUAAIBBxBUAAIBBxBUAAIBBxBUAAIBBxBUAAIBBc/KHmwHMrMW7j0/6vm/vXW1wJgAw+7FyBQAAYBBxBQAAYBBxBQAAYBBxBQAAYBBxBQAAYBBxBQAAYBBxBQAAYBBxBQAAYBBxBQAAYBBxBQAAYBBxBQAAYBBxBQAAYFDccfXHP/5RX/3qV5Wbm6v09HR95jOfUXd3t33csizV19fL6/UqLS1NZWVl6unpiXmMSCSimpoa5eXlKSMjQ2vXrlV/f//UzwYAACDB4oqroaEh3X333UpJSdFzzz2n1157TU8++aRuvvlme0xjY6OamprU3Nysrq4ueTweVVZWanh42B7j9/vV1tam1tZWnTp1SiMjI1qzZo3GxsaMnRgAAEAiOCzLsj7q4N27d+tXv/qVXn755SsetyxLXq9Xfr9fjzzyiKT3V6ncbrf27dunzZs3KxQKad68eTpy5IjWr18vSbpw4YJ8Pp9OnDihlStXfug8wuGwXC6XQqGQsrKyPur0AUzS4t3HJ33ft/euNjgTAEiMeNojrpWrZ599VsuWLdOXv/xlzZ8/X0uXLtVTTz1lH+/t7VUwGFRVVZW9z+l0asWKFers7JQkdXd369KlSzFjvF6vCgsL7THjRSIRhcPhmA0AAGA2iiuu3nrrLR08eFAFBQV64YUXtGXLFn3rW9/SM888I0kKBoOSJLfbHXM/t9ttHwsGg0pNTVV2dvZVx4zX0NAgl8tlbz6fL55pAwAAzJi44ury5cu68847FQgEtHTpUm3evFnf+MY3dPDgwZhxDocj5rZlWRP2jXetMXV1dQqFQvbW19cXz7QBAABmTFxxtWDBAt12220x+2699VadP39ekuTxeCRpwgrU4OCgvZrl8XgUjUY1NDR01THjOZ1OZWVlxWwAAACzUVxxdffdd+v111+P2ff73/9eixYtkiTl5+fL4/Govb3dPh6NRtXR0aHS0lJJUnFxsVJSUmLGDAwM6Ny5c/YYAACAuSo5nsEPP/ywSktLFQgEVF1drV//+tc6dOiQDh06JOn9twP9fr8CgYAKCgpUUFCgQCCg9PR0bdiwQZLkcrm0adMm7dy5U7m5ucrJyVFtba2KiopUUVFh/gwBAABmUFxxddddd6mtrU11dXV6/PHHlZ+frwMHDujBBx+0x+zatUujo6PaunWrhoaGVFJSopMnTyozM9Mes3//fiUnJ6u6ulqjo6MqLy9XS0uLkpKSzJ0ZgBhTuZwCAOCji+s6V7MF17kC4peouOI6VwCuB9N2nSsAAABcG3EFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgUFxxVV9fL4fDEbN5PB77uGVZqq+vl9frVVpamsrKytTT0xPzGJFIRDU1NcrLy1NGRobWrl2r/v5+M2cDAACQYHGvXN1+++0aGBiwt7Nnz9rHGhsb1dTUpObmZnV1dcnj8aiyslLDw8P2GL/fr7a2NrW2turUqVMaGRnRmjVrNDY2ZuaMAAAAEig57jskJ8esVn3AsiwdOHBAe/bs0bp16yRJhw8fltvt1rFjx7R582aFQiE9/fTTOnLkiCoqKiRJR48elc/n04svvqiVK1dO8XQAAAASK+6VqzfeeENer1f5+fl64IEH9NZbb0mSent7FQwGVVVVZY91Op1asWKFOjs7JUnd3d26dOlSzBiv16vCwkJ7zJVEIhGFw+GYDQAAYDaKK65KSkr0zDPP6IUXXtBTTz2lYDCo0tJSvfvuuwoGg5Ikt9sdcx+3220fCwaDSk1NVXZ29lXHXElDQ4NcLpe9+Xy+eKYNAAAwY+KKq1WrVun+++9XUVGRKioqdPz4cUnvv/33AYfDEXMfy7Im7Bvvw8bU1dUpFArZW19fXzzTBgAAmDFTuhRDRkaGioqK9MYbb9ifwxq/AjU4OGivZnk8HkWjUQ0NDV11zJU4nU5lZWXFbAAAALNR3B9o/78ikYh+97vf6Z577lF+fr48Ho/a29u1dOlSSVI0GlVHR4f27dsnSSouLlZKSora29tVXV0tSRoYGNC5c+fU2Ng4xVMBrn+Ldx9P9BQAAB8irriqra3Vvffeq4ULF2pwcFBPPPGEwuGwNm7cKIfDIb/fr0AgoIKCAhUUFCgQCCg9PV0bNmyQJLlcLm3atEk7d+5Ubm6ucnJyVFtba7/NCAAAMNfFFVf9/f36yle+onfeeUfz5s3TZz/7WZ0+fVqLFi2SJO3atUujo6PaunWrhoaGVFJSopMnTyozM9N+jP379ys5OVnV1dUaHR1VeXm5WlpalJSUZPbMAAAAEsBhWZaV6EnEKxwOy+VyKRQK8fkr3FDm4tuCb+9dnegpAMCUxdMe/LYgAACAQcQVAACAQcQVAACAQVO6FAMAfJipfE6Mz2sBmItYuQIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADCIuAIAADBoSnHV0NAgh8Mhv99v77MsS/X19fJ6vUpLS1NZWZl6enpi7heJRFRTU6O8vDxlZGRo7dq16u/vn8pUAAAAZoVJx1VXV5cOHTqkJUuWxOxvbGxUU1OTmpub1dXVJY/Ho8rKSg0PD9tj/H6/2tra1NraqlOnTmlkZERr1qzR2NjY5M8EAABgFphUXI2MjOjBBx/UU089pezsbHu/ZVk6cOCA9uzZo3Xr1qmwsFCHDx/We++9p2PHjkmSQqGQnn76aT355JOqqKjQ0qVLdfToUZ09e1YvvvjiFf9eJBJROByO2QAAAGajScXVtm3btHr1alVUVMTs7+3tVTAYVFVVlb3P6XRqxYoV6uzslCR1d3fr0qVLMWO8Xq8KCwvtMeM1NDTI5XLZm8/nm8y0AQAApl3ccdXa2qrf/OY3amhomHAsGAxKktxud8x+t9ttHwsGg0pNTY1Z8Ro/Zry6ujqFQiF76+vri3faAAAAMyI5nsF9fX369re/rZMnT+qmm2666jiHwxFz27KsCfvGu9YYp9Mpp9MZz1QBAAASIq6Vq+7ubg0ODqq4uFjJyclKTk5WR0eHvv/97ys5OdlesRq/AjU4OGgf83g8ikajGhoauuoYAACAuSquuCovL9fZs2d15swZe1u2bJkefPBBnTlzRh//+Mfl8XjU3t5u3ycajaqjo0OlpaWSpOLiYqWkpMSMGRgY0Llz5+wxAAAAc1VcbwtmZmaqsLAwZl9GRoZyc3Pt/X6/X4FAQAUFBSooKFAgEFB6ero2bNggSXK5XNq0aZN27typ3Nxc5eTkqLa2VkVFRRM+IA8AADDXxBVXH8WuXbs0OjqqrVu3amhoSCUlJTp58qQyMzPtMfv371dycrKqq6s1Ojqq8vJytbS0KCkpyfR0AAAAZpTDsiwr0ZOIVzgclsvlUigUUlZWVqKnA8yYxbuPJ3oKM+rtvasTPQUAkBRfe/DbggAAAAYZf1sQAEyZ6kodK18AEoG4AmbYjfbWHgDcaHhbEAAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwKC44urgwYNasmSJsrKylJWVpeXLl+u5556zj1uWpfr6enm9XqWlpamsrEw9PT0xjxGJRFRTU6O8vDxlZGRo7dq16u/vN3M2AAAACRZXXN1yyy3au3evXnnlFb3yyiv64he/qC996Ut2QDU2NqqpqUnNzc3q6uqSx+NRZWWlhoeH7cfw+/1qa2tTa2urTp06pZGREa1Zs0ZjY2NmzwwAACABHJZlWVN5gJycHH33u9/V17/+dXm9Xvn9fj3yyCOS3l+lcrvd2rdvnzZv3qxQKKR58+bpyJEjWr9+vSTpwoUL8vl8OnHihFauXPmR/mY4HJbL5VIoFFJWVtZUpg/MuMW7jyd6CjeMt/euTvQUAFwn4mmPSX/mamxsTK2trbp48aKWL1+u3t5eBYNBVVVV2WOcTqdWrFihzs5OSVJ3d7cuXboUM8br9aqwsNAecyWRSEThcDhmAwAAmI3ijquzZ8/qYx/7mJxOp7Zs2aK2tjbddtttCgaDkiS32x0z3u1228eCwaBSU1OVnZ191TFX0tDQIJfLZW8+ny/eaQMAAMyIuOPqU5/6lM6cOaPTp0/rm9/8pjZu3KjXXnvNPu5wOGLGW5Y1Yd94Hzamrq5OoVDI3vr6+uKdNgAAwIyIO65SU1P1iU98QsuWLVNDQ4PuuOMOfe9735PH45GkCStQg4OD9mqWx+NRNBrV0NDQVcdcidPptL+h+MEGAAAwG035OleWZSkSiSg/P18ej0ft7e32sWg0qo6ODpWWlkqSiouLlZKSEjNmYGBA586ds8cAAADMZcnxDP7Od76jVatWyefzaXh4WK2trfrlL3+p559/Xg6HQ36/X4FAQAUFBSooKFAgEFB6ero2bNggSXK5XNq0aZN27typ3Nxc5eTkqLa2VkVFRaqoqJiWEwQAAJhJccXVn/70Jz300EMaGBiQy+XSkiVL9Pzzz6uyslKStGvXLo2Ojmrr1q0aGhpSSUmJTp48qczMTPsx9u/fr+TkZFVXV2t0dFTl5eVqaWlRUlKS2TMDAABIgClf5yoRuM4V5jKuczVzuM4VAFNm5DpXAAAAmIi4AgAAMIi4AgAAMIi4AgAAMIi4AgAAMIi4AgAAMIi4AgAAMIi4AgAAMCiuK7QDeB8XAgUAXA0rVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYlJ3oCADBdFu8+Pun7vr13tcGZALiRsHIFAABgUFxx1dDQoLvuukuZmZmaP3++7rvvPr3++usxYyzLUn19vbxer9LS0lRWVqaenp6YMZFIRDU1NcrLy1NGRobWrl2r/v7+qZ8NAABAgsUVVx0dHdq2bZtOnz6t9vZ2/fWvf1VVVZUuXrxoj2lsbFRTU5Oam5vV1dUlj8ejyspKDQ8P22P8fr/a2trU2tqqU6dOaWRkRGvWrNHY2Ji5MwMAAEgAh2VZ1mTv/Oc//1nz589XR0eHPv/5z8uyLHm9Xvn9fj3yyCOS3l+lcrvd2rdvnzZv3qxQKKR58+bpyJEjWr9+vSTpwoUL8vl8OnHihFauXPmhfzccDsvlcikUCikrK2uy0wcmbSqf5cHcwGeuAPxf8bTHlD5zFQqFJEk5OTmSpN7eXgWDQVVVVdljnE6nVqxYoc7OTklSd3e3Ll26FDPG6/WqsLDQHjNeJBJROByO2QAAAGajSceVZVnasWOHPve5z6mwsFCSFAwGJUlutztmrNvtto8Fg0GlpqYqOzv7qmPGa2hokMvlsjefzzfZaQMAAEyrScfV9u3b9dvf/lY/+tGPJhxzOBwxty3LmrBvvGuNqaurUygUsre+vr7JThsAAGBaTSquampq9Oyzz+oXv/iFbrnlFnu/x+ORpAkrUIODg/ZqlsfjUTQa1dDQ0FXHjOd0OpWVlRWzAQAAzEZxxZVlWdq+fbt+8pOf6Oc//7ny8/Njjufn58vj8ai9vd3eF41G1dHRodLSUklScXGxUlJSYsYMDAzo3Llz9hgAAIC5Kq4rtG/btk3Hjh3Tz372M2VmZtorVC6XS2lpaXI4HPL7/QoEAiooKFBBQYECgYDS09O1YcMGe+ymTZu0c+dO5ebmKicnR7W1tSoqKlJFRYX5MwQAAJhBccXVwYMHJUllZWUx+3/4wx/qa1/7miRp165dGh0d1datWzU0NKSSkhKdPHlSmZmZ9vj9+/crOTlZ1dXVGh0dVXl5uVpaWpSUlDS1swEAAEiwKV3nKlG4zhUSjetcXf+4zhWA/2vGrnMFAACAWMQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQXH9/A1wPeEq6wCA6cDKFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHJiZ4AAMxGi3cfn/R939672uBMAMw1ca9cvfTSS7r33nvl9XrlcDj005/+NOa4ZVmqr6+X1+tVWlqaysrK1NPTEzMmEomopqZGeXl5ysjI0Nq1a9Xf3z+lEwEAAJgN4o6rixcv6o477lBzc/MVjzc2NqqpqUnNzc3q6uqSx+NRZWWlhoeH7TF+v19tbW1qbW3VqVOnNDIyojVr1mhsbGzyZwIAADALxP224KpVq7Rq1aorHrMsSwcOHNCePXu0bt06SdLhw4fldrt17Ngxbd68WaFQSE8//bSOHDmiiooKSdLRo0fl8/n04osvauXKlVM4HQAAgMQy+oH23t5eBYNBVVVV2fucTqdWrFihzs5OSVJ3d7cuXboUM8br9aqwsNAeM14kElE4HI7ZAAAAZiOjcRUMBiVJbrc7Zr/b7baPBYNBpaamKjs7+6pjxmtoaJDL5bI3n89nctoAAADGTMu3BR0OR8xty7Im7BvvWmPq6uq0Y8cO+3Y4HCawMKVvcwEAMF2Mrlx5PB5JmrACNTg4aK9meTweRaNRDQ0NXXXMeE6nU1lZWTEbAADAbGQ0rvLz8+XxeNTe3m7vi0aj6ujoUGlpqSSpuLhYKSkpMWMGBgZ07tw5ewwAAMBcFffbgiMjI3rzzTft2729vTpz5oxycnK0cOFC+f1+BQIBFRQUqKCgQIFAQOnp6dqwYYMkyeVyadOmTdq5c6dyc3OVk5Oj2tpaFRUV2d8eBAAAmKvijqtXXnlFX/jCF+zbH3wWauPGjWppadGuXbs0OjqqrVu3amhoSCUlJTp58qQyMzPt++zfv1/Jycmqrq7W6OioysvL1dLSoqSkJAOnBAAAkDgOy7KsRE8iXuFwWC6XS6FQiM9f3cD4QDtmK37+Brj+xNMe/HAzAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQXH/tiAA4Nqm8tNM/HQOMPexcgUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQP3+DhJrKz4QAADAbsXIFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEHEFAABgEJdiAIBZZCqXJ3l772qDMwEwWcQVpoxrVQEA8L94WxAAAMAg4goAAMAg4goAAMAg4goAAMAg4goAAMAgvi0ISXzjD7gecBkHYHZI6MrVD37wA+Xn5+umm25ScXGxXn755UROBwAAYMoSFlc//vGP5ff7tWfPHr366qu65557tGrVKp0/fz5RUwIAAJgyh2VZViL+cElJie68804dPHjQ3nfrrbfqvvvuU0NDwzXvGw6H5XK5FAqFlJWVNd1TnTN4aw/AZPG2IHBt8bRHQj5zFY1G1d3drd27d8fsr6qqUmdn54TxkUhEkUjEvh0KhSS9f6L4X5cj7yV6CgDmqIUP/7+E/e1zj61M2N8GPqoPmuOjrEklJK7eeecdjY2Nye12x+x3u90KBoMTxjc0NOixxx6bsN/n803bHAEAM8N1INEzAD664eFhuVyua45J6LcFHQ5HzG3Lsibsk6S6ujrt2LHDvn358mX95S9/UW5u7hXHY3LC4bB8Pp/6+vp4u3Ua8TxPP57jmcHzPDN4nmfGhz3PlmVpeHhYXq/3Qx8rIXGVl5enpKSkCatUg4ODE1azJMnpdMrpdMbsu/nmm6dzije0rKws/gWeATzP04/neGbwPM8MnueZca3n+cNWrD6QkG8Lpqamqri4WO3t7TH729vbVVpamogpAQAAGJGwtwV37Nihhx56SMuWLdPy5ct16NAhnT9/Xlu2bEnUlAAAAKYsYXG1fv16vfvuu3r88cc1MDCgwsJCnThxQosWLUrUlG54TqdTjz766IS3YGEWz/P04zmeGTzPM4PneWaYfJ4Tdp0rAACA6xE/3AwAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGAQcQUAAGBQcqInMBmXL1/WhQsXlJmZKYfDkejpAACA65xlWRoeHpbX69Xf/M2116bmZFxduHBBPp8v0dMAAAA3mL6+Pt1yyy3XHDMn4yozM1PS+yeYlZWV4NkAAIDrXTgcls/nsxvkWuZkXH3wVmBWVhZxBQAAZsxH+TgSH2gHAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwaE7+cDOAxFi8+/i0/423966e9r8BANOJlSsAAACDiCsAAACDiCsAAACD+MwVgFmFz3UBmOtYuQIAADCIuAIAADCIuAIAADCIuAIAADAo7rj64x//qK9+9avKzc1Venq6PvOZz6i7u9s+blmW6uvr5fV6lZaWprKyMvX09MQ8RiQSUU1NjfLy8pSRkaG1a9eqv79/6mcDAACQYHHF1dDQkO6++26lpKToueee02uvvaYnn3xSN998sz2msbFRTU1Nam5uVldXlzwejyorKzU8PGyP8fv9amtrU2trq06dOqWRkRGtWbNGY2Njxk4MAAAgERyWZVkfdfDu3bv1q1/9Si+//PIVj1uWJa/XK7/fr0ceeUTS+6tUbrdb+/bt0+bNmxUKhTRv3jwdOXJE69evlyRduHBBPp9PJ06c0MqVKz90HuFwWC6XS6FQSFlZWR91+gCmaCYukzATuBQDgHjF0x5xrVw9++yzWrZsmb785S9r/vz5Wrp0qZ566in7eG9vr4LBoKqqqux9TqdTK1asUGdnpySpu7tbly5dihnj9XpVWFhojxkvEokoHA7HbAAAALNRXHH11ltv6eDBgyooKNALL7ygLVu26Fvf+paeeeYZSVIwGJQkud3umPu53W77WDAYVGpqqrKzs686ZryGhga5XC578/l88UwbAABgxsQVV5cvX9add96pQCCgpUuXavPmzfrGN76hgwcPxoxzOBwxty3LmrBvvGuNqaurUygUsre+vr54pg0AADBj4oqrBQsW6LbbbovZd+utt+r8+fOSJI/HI0kTVqAGBwft1SyPx6NoNKqhoaGrjhnP6XQqKysrZgMAAJiN4oqru+++W6+//nrMvt///vdatGiRJCk/P18ej0ft7e328Wg0qo6ODpWWlkqSiouLlZKSEjNmYGBA586ds8cAAADMVXH9cPPDDz+s0tJSBQIBVVdX69e//rUOHTqkQ4cOSXr/7UC/369AIKCCggIVFBQoEAgoPT1dGzZskCS5XC5t2rRJO3fuVG5urnJyclRbW6uioiJVVFSYP0MAAIAZFFdc3XXXXWpra1NdXZ0ef/xx5efn68CBA3rwwQftMbt27dLo6Ki2bt2qoaEhlZSU6OTJk8rMzLTH7N+/X8nJyaqurtbo6KjKy8vV0tKipKQkc2cGAACQAHFd52q24DpXQGJwnSsAN6ppu84VAAAAro24AgAAMIi4AgAAMIi4AgAAMIi4AgAAMCiuSzEAmN2ul2/zAcBcxsoVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQXHFVX19vRwOR8zm8Xjs45Zlqb6+Xl6vV2lpaSorK1NPT0/MY0QiEdXU1CgvL08ZGRlau3at+vv7zZwNAABAgsW9cnX77bdrYGDA3s6ePWsfa2xsVFNTk5qbm9XV1SWPx6PKykoNDw/bY/x+v9ra2tTa2qpTp05pZGREa9as0djYmJkzAgAASKDkuO+QnByzWvUBy7J04MAB7dmzR+vWrZMkHT58WG63W8eOHdPmzZsVCoX09NNP68iRI6qoqJAkHT16VD6fTy+++KJWrlw5xdMBAABIrLhXrt544w15vV7l5+frgQce0FtvvSVJ6u3tVTAYVFVVlT3W6XRqxYoV6uzslCR1d3fr0qVLMWO8Xq8KCwvtMVcSiUQUDodjNgAAgNkorrgqKSnRM888oxdeeEFPPfWUgsGgSktL9e677yoYDEqS3G53zH3cbrd9LBgMKjU1VdnZ2VcdcyUNDQ1yuVz25vP54pk2AADAjIkrrlatWqX7779fRUVFqqio0PHjxyW9//bfBxwOR8x9LMuasG+8DxtTV1enUChkb319ffFMGwAAYMZM6VIMGRkZKioq0htvvGF/Dmv8CtTg4KC9muXxeBSNRjU0NHTVMVfidDqVlZUVswEAAMxGU4qrSCSi3/3ud1qwYIHy8/Pl8XjU3t5uH49Go+ro6FBpaakkqbi4WCkpKTFjBgYGdO7cOXsMAADAXBbXtwVra2t17733auHChRocHNQTTzyhcDisjRs3yuFwyO/3KxAIqKCgQAUFBQoEAkpPT9eGDRskSS6XS5s2bdLOnTuVm5urnJwc1dbW2m8zAgAAzHVxxVV/f7++8pWv6J133tG8efP02c9+VqdPn9aiRYskSbt27dLo6Ki2bt2qoaEhlZSU6OTJk8rMzLQfY//+/UpOTlZ1dbVGR0dVXl6ulpYWJSUlmT0zAACABHBYlmUlehLxCofDcrlcCoVCfP4K+D8W7z6e6CnMCW/vXZ3oKQCYY+JpD35bEAAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwCDiCgAAwKDkRE8AuFEs3n080VMAAMwAVq4AAAAMIq4AAAAMIq4AAAAMIq4AAAAMIq4AAAAMIq4AAAAMIq4AAAAMIq4AAAAMIq4AAAAMIq4AAAAM4udvANxwpvuniN7eu3paHx/A7MbKFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEFTiquGhgY5HA75/X57n2VZqq+vl9frVVpamsrKytTT0xNzv0gkopqaGuXl5SkjI0Nr165Vf3//VKYCAAAwK0w6rrq6unTo0CEtWbIkZn9jY6OamprU3Nysrq4ueTweVVZWanh42B7j9/vV1tam1tZWnTp1SiMjI1qzZo3GxsYmfyYAAACzwKTiamRkRA8++KCeeuopZWdn2/sty9KBAwe0Z88erVu3ToWFhTp8+LDee+89HTt2TJIUCoX09NNP68knn1RFRYWWLl2qo0eP6uzZs3rxxRfNnBUAAECCTCqutm3bptWrV6uioiJmf29vr4LBoKqqqux9TqdTK1asUGdnpySpu7tbly5dihnj9XpVWFhojxkvEokoHA7HbAAAALNRcrx3aG1t1W9+8xt1dXVNOBYMBiVJbrc7Zr/b7dYf/vAHe0xqamrMitcHYz64/3gNDQ167LHH4p0qAADAjItr5aqvr0/f/va3dfToUd10001XHedwOGJuW5Y1Yd941xpTV1enUChkb319ffFMGwAAYMbEFVfd3d0aHBxUcXGxkpOTlZycrI6ODn3/+99XcnKyvWI1fgVqcHDQPubxeBSNRjU0NHTVMeM5nU5lZWXFbAAAALNRXHFVXl6us2fP6syZM/a2bNkyPfjggzpz5ow+/vGPy+PxqL293b5PNBpVR0eHSktLJUnFxcVKSUmJGTMwMKBz587ZYwAAAOaquD5zlZmZqcLCwph9GRkZys3Ntff7/X4FAgEVFBSooKBAgUBA6enp2rBhgyTJ5XJp06ZN2rlzp3Jzc5WTk6Pa2loVFRVN+IA8AADAXBP3B9o/zK5duzQ6OqqtW7dqaGhIJSUlOnnypDIzM+0x+/fvV3JysqqrqzU6Oqry8nK1tLQoKSnJ9HQAAABmlMOyLCvRk4hXOByWy+VSKBTi81eYMxbvPp7oKWCGvL13daKnAMCweNqD3xYEAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwiLgCAAAwKK64OnjwoJYsWaKsrCxlZWVp+fLleu655+zjlmWpvr5eXq9XaWlpKisrU09PT8xjRCIR1dTUKC8vTxkZGVq7dq36+/vNnA0AAECCxRVXt9xyi/bu3atXXnlFr7zyir74xS/qS1/6kh1QjY2NampqUnNzs7q6uuTxeFRZWanh4WH7Mfx+v9ra2tTa2qpTp05pZGREa9as0djYmNkzAwAASACHZVnWVB4gJydH3/3ud/X1r39dXq9Xfr9fjzzyiKT3V6ncbrf27dunzZs3KxQKad68eTpy5IjWr18vSbpw4YJ8Pp9OnDihlStXfqS/GQ6H5XK5FAqFlJWVNZXpAzNm8e7jiZ4CZsjbe1cnegoADIunPSb9mauxsTG1trbq4sWLWr58uXp7exUMBlVVVWWPcTqdWrFihTo7OyVJ3d3dunTpUswYr9erwsJCe8yVRCIRhcPhmA0AAGA2ijuuzp49q4997GNyOp3asmWL2tradNtttykYDEqS3G53zHi3220fCwaDSk1NVXZ29lXHXElDQ4NcLpe9+Xy+eKcNAAAwI+KOq0996lM6c+aMTp8+rW9+85vauHGjXnvtNfu4w+GIGW9Z1oR9433YmLq6OoVCIXvr6+uLd9oAAAAzIjneO6SmpuoTn/iEJGnZsmXq6urS9773PftzVsFgUAsWLLDHDw4O2qtZHo9H0WhUQ0NDMatXg4ODKi0tverfdDqdcjqd8U4VABJiJj5fx+e6gNlryte5sixLkUhE+fn58ng8am9vt49Fo1F1dHTY4VRcXKyUlJSYMQMDAzp37tw14woAAGCuiGvl6jvf+Y5WrVoln8+n4eFhtba26pe//KWef/55ORwO+f1+BQIBFRQUqKCgQIFAQOnp6dqwYYMkyeVyadOmTdq5c6dyc3OVk5Oj2tpaFRUVqaKiYlpOEAAAYCbFFVd/+tOf9NBDD2lgYEAul0tLlizR888/r8rKSknSrl27NDo6qq1bt2poaEglJSU6efKkMjMz7cfYv3+/kpOTVV1drdHRUZWXl6ulpUVJSUlmzwwAACABpnydq0TgOleYi7jOFUziM1fAzJqR61wBAABgIuIKAADAIOIKAADAIOIKAADAIOIKAADAIOIKAADAIOIKAADAIOIKAADAIOIKAADAIOIKAADAIOIKAADAIOIKAADAIOIKAADAIOIKAADAIOIKAADAIOIKAADAIOIKAADAIOIKAADAIOIKAADAIOIKAADAIOIKAADAIOIKAADAIOIKAADAoORETwCYDRbvPp7oKQAArhOsXAEAABhEXAEAABhEXAEAABhEXAEAABjEB9oBYA6aiS9hvL139bT/DeB6xMoVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQcQVAACAQXHFVUNDg+666y5lZmZq/vz5uu+++/T666/HjLEsS/X19fJ6vUpLS1NZWZl6enpixkQiEdXU1CgvL08ZGRlau3at+vv7p342AAAACRZXXHV0dGjbtm06ffq02tvb9de//lVVVVW6ePGiPaaxsVFNTU1qbm5WV1eXPB6PKisrNTw8bI/x+/1qa2tTa2urTp06pZGREa1Zs0ZjY2PmzgwAACABHJZlWZO985///GfNnz9fHR0d+vznPy/LsuT1euX3+/XII49Ien+Vyu12a9++fdq8ebNCoZDmzZunI0eOaP369ZKkCxcuyOfz6cSJE1q5cuWH/t1wOCyXy6VQKKSsrKzJTh+wzcRPiQBzDT9/A/yveNpjSp+5CoVCkqScnBxJUm9vr4LBoKqqquwxTqdTK1asUGdnpySpu7tbly5dihnj9XpVWFhojxkvEokoHA7HbAAAALPRpOPKsizt2LFDn/vc51RYWChJCgaDkiS32x0z1u1228eCwaBSU1OVnZ191THjNTQ0yOVy2ZvP55vstAEAAKbVpONq+/bt+u1vf6sf/ehHE445HI6Y25ZlTdg33rXG1NXVKRQK2VtfX99kpw0AADCtJhVXNTU1evbZZ/WLX/xCt9xyi73f4/FI0oQVqMHBQXs1y+PxKBqNamho6KpjxnM6ncrKyorZAAAAZqO44sqyLG3fvl0/+clP9POf/1z5+fkxx/Pz8+XxeNTe3m7vi0aj6ujoUGlpqSSpuLhYKSkpMWMGBgZ07tw5ewwAAMBclRzP4G3btunYsWP62c9+pszMTHuFyuVyKS0tTQ6HQ36/X4FAQAUFBSooKFAgEFB6ero2bNhgj920aZN27typ3Nxc5eTkqLa2VkVFRaqoqDB/hgAAADMorrg6ePCgJKmsrCxm/w9/+EN97WtfkyTt2rVLo6Oj2rp1q4aGhlRSUqKTJ08qMzPTHr9//34lJyerurpao6OjKi8vV0tLi5KSkqZ2NgAAAAk2petcJQrXuYJpXOcKmIjrXAH/a8aucwUAAIBYxBUAAIBBxBUAAIBBxBUAAIBBxBUAAIBBxBUAAIBBxBUAAIBBxBUAAIBBxBUAAIBBxBUAAIBBxBUAAIBBxBUAAIBBxBUAAIBBxBUAAIBBxBUAAIBBxBUAAIBBxBUAAIBByYmeAABgdlq8+/i0Pv7be1dP6+MDicLKFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHEFQAAgEHJiZ4AAODGtHj38Wn/G2/vXT3tfwMYL+6Vq5deekn33nuvvF6vHA6HfvrTn8YctyxL9fX18nq9SktLU1lZmXp6emLGRCIR1dTUKC8vTxkZGVq7dq36+/undCIAAACzQdwrVxcvXtQdd9yhf/qnf9L9998/4XhjY6OamprU0tKiT37yk3riiSdUWVmp119/XZmZmZIkv9+v//qv/1Jra6tyc3O1c+dOrVmzRt3d3UpKSpr6WeG6MxP/hwsAgAlxx9WqVau0atWqKx6zLEsHDhzQnj17tG7dOknS4cOH5Xa7dezYMW3evFmhUEhPP/20jhw5ooqKCknS0aNH5fP59OKLL2rlypUTHjcSiSgSidi3w+FwvNMGAACYEUY/0N7b26tgMKiqqip7n9Pp1IoVK9TZ2SlJ6u7u1qVLl2LGeL1eFRYW2mPGa2hokMvlsjefz2dy2gAAAMYYjatgMChJcrvdMfvdbrd9LBgMKjU1VdnZ2VcdM15dXZ1CoZC99fX1mZw2AACAMdPybUGHwxFz27KsCfvGu9YYp9Mpp9NpbH4AAADTxejKlcfjkaQJK1CDg4P2apbH41E0GtXQ0NBVxwAAAMxVRuMqPz9fHo9H7e3t9r5oNKqOjg6VlpZKkoqLi5WSkhIzZmBgQOfOnbPHAAAAzFVxvy04MjKiN998077d29urM2fOKCcnRwsXLpTf71cgEFBBQYEKCgoUCASUnp6uDRs2SJJcLpc2bdqknTt3Kjc3Vzk5OaqtrVVRUZH97UEAAIC5Ku64euWVV/SFL3zBvr1jxw5J0saNG9XS0qJdu3ZpdHRUW7du1dDQkEpKSnTy5En7GleStH//fiUnJ6u6ulqjo6MqLy9XS0sL17gCAABznsOyLCvRk4hXOByWy+VSKBRSVlZWoqeDGcBFRAFMBj9/A1PiaQ9+uBkAAMAg4goAAMAg4goAAMAg4goAAMAg4goAAMAg4goAAMAg4goAAMAg4goAAMCguK/QDgDAXDETFyDmQqUYj5UrAAAAg4grAAAAg4grAAAAg4grAAAAg4grAAAAg4grAAAAg4grAAAAg4grAAAAg4grAAAAg4grAAAAg4grAAAAg/htQUzZTPx2FwAAcwUrVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAZxnSsAAKZguq/19/be1dP6+DCPlSsAAACDiCsAAACDiCsAAACDiCsAAACDiCsAAACDiCsAAACDuBTDdW66vyIMAABiJTSufvCDH+i73/2uBgYGdPvtt+vAgQO65557EjklAABmlZn4n2SupWVWwt4W/PGPfyy/3689e/bo1Vdf1T333KNVq1bp/PnziZoSAADAlCUsrpqamrRp0yb98z//s2699VYdOHBAPp9PBw8eTNSUAAAApiwhbwtGo1F1d3dr9+7dMfurqqrU2dk5YXwkElEkErFvh0IhSVI4HJ7eic6AwkdfSPQUAAA3uOvhv6fT7YPnyLKsDx2bkLh65513NDY2JrfbHbPf7XYrGAxOGN/Q0KDHHntswn6fzzdtcwQA4EbhOpDoGcwdw8PDcrlc1xyT0A+0OxyOmNuWZU3YJ0l1dXXasWOHffvy5cv6y1/+otzc3CuOx+SEw2H5fD719fUpKysr0dO5bvE8Tz+e45nB8zwzeJ5nxoc9z5ZlaXh4WF6v90MfKyFxlZeXp6SkpAmrVIODgxNWsyTJ6XTK6XTG7Lv55punc4o3tKysLP4FngE8z9OP53hm8DzPDJ7nmXGt5/nDVqw+kJAPtKempqq4uFjt7e0x+9vb21VaWpqIKQEAABiRsLcFd+zYoYceekjLli3T8uXLdejQIZ0/f15btmxJ1JQAAACmLGFxtX79er377rt6/PHHNTAwoMLCQp04cUKLFi1K1JRueE6nU48++uiEt2BhFs/z9OM5nhk8zzOD53lmmHyeHdZH+U4hAAAAPhJ+uBkAAMAg4goAAMAg4goAAMAg4goAAMAg4goAAMAg4gpXtHjxYjkcjpht/A9tI34/+MEPlJ+fr5tuuknFxcV6+eWXEz2l60p9ff2E163H40n0tOa8l156Sffee6+8Xq8cDod++tOfxhy3LEv19fXyer1KS0tTWVmZenp6EjPZOezDnuevfe1rE17fn/3sZxMz2TmqoaFBd911lzIzMzV//nzdd999ev3112PGmHg9E1e4qg+uQfbB9i//8i+JntKc9uMf/1h+v1979uzRq6++qnvuuUerVq3S+fPnEz2168rtt98e87o9e/Zsoqc05128eFF33HGHmpubr3i8sbFRTU1Nam5uVldXlzwejyorKzU8PDzDM53bPux5lqS///u/j3l9nzhxYgZnOPd1dHRo27ZtOn36tNrb2/XXv/5VVVVVunjxoj3GyOvZAq5g0aJF1v79+xM9jevK3/3d31lbtmyJ2ffpT3/a2r17d4JmdP159NFHrTvuuCPR07iuSbLa2trs25cvX7Y8Ho+1d+9ee9///M//WC6Xy/qP//iPBMzw+jD+ebYsy9q4caP1pS99KSHzuV4NDg5akqyOjg7Lssy9nlm5wlXt27dPubm5+sxnPqN/+7d/UzQaTfSU5qxoNKru7m5VVVXF7K+qqlJnZ2eCZnV9euONN+T1epWfn68HHnhAb731VqKndF3r7e1VMBiMeW07nU6tWLGC1/Y0+OUvf6n58+frk5/8pL7xjW9ocHAw0VOa00KhkCQpJydHkrnXc8J+/gaz27e//W3deeedys7O1q9//WvV1dWpt7dX//mf/5noqc1J77zzjsbGxuR2u2P2u91uBYPBBM3q+lNSUqJnnnlGn/zkJ/WnP/1JTzzxhEpLS9XT06Pc3NxET++69MHr90qv7T/84Q+JmNJ1a9WqVfryl7+sRYsWqbe3V//6r/+qL37xi+ru7uancSbBsizt2LFDn/vc51RYWCjJ3OuZuLqB1NfX67HHHrvmmK6uLi1btkwPP/ywvW/JkiXKzs7WP/7jP9qrWZgch8MRc9uyrAn7MHmrVq2y/7moqEjLly/X3/7t3+rw4cPasWNHAmd2/eO1Pf3Wr19v/3NhYaGWLVumRYsW6fjx41q3bl0CZzY3bd++Xb/97W916tSpCcem+nomrm4g27dv1wMPPHDNMYsXL77i/g++kfLmm28SV5OQl5enpKSkCatUg4ODE/4PCeZkZGSoqKhIb7zxRqKnct364NuYwWBQCxYssPfz2p5+CxYs0KJFi3h9T0JNTY2effZZvfTSS7rlllvs/aZez3zm6gaSl5enT3/609fcbrrppive99VXX5WkmBcbPrrU1FQVFxervb09Zn97e7tKS0sTNKvrXyQS0e9+9ztet9MoPz9fHo8n5rUdjUbV0dHBa3uavfvuu+rr6+P1HQfLsrR9+3b95Cc/0c9//nPl5+fHHDf1emblChP893//t06fPq0vfOELcrlc6urq0sMPP6y1a9dq4cKFiZ7enLVjxw499NBDWrZsmZYvX65Dhw7p/Pnz2rJlS6Kndt2ora3Vvffeq4ULF2pwcFBPPPGEwuGwNm7cmOipzWkjIyN688037du9vb06c+aMcnJytHDhQvn9fgUCARUUFKigoECBQEDp6enasGFDAmc991zrec7JyVF9fb3uv/9+LViwQG+//ba+853vKC8vT//wD/+QwFnPLdu2bdOxY8f0s5/9TJmZmfa7CS6XS2lpaXI4HGZezya/0ojrQ3d3t1VSUmK5XC7rpptusj71qU9Zjz76qHXx4sVET23O+/d//3dr0aJFVmpqqnXnnXfaX/+FGevXr7cWLFhgpaSkWF6v11q3bp3V09OT6GnNeb/4xS8sSRO2jRs3Wpb1/tfXH330Ucvj8VhOp9P6/Oc/b509ezaxk56DrvU8v/fee1ZVVZU1b948KyUlxVq4cKG1ceNG6/z584me9pxypedXkvXDH/7QHmPi9ez4//8YAAAADOAzVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAYRVwAAAAb9f2CaSAAv5NHyAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 700x900 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(2, 1, figsize=(7, 9), sharex=True)\n",
    "\n",
    "y1_grs = p1.posterior[\"mu\"].sel(__obs__=0).to_numpy().flatten()\n",
    "y2_grs = p1.posterior[\"mu\"].sel(__obs__=1).to_numpy().flatten()\n",
    "\n",
    "axes[0].hist(y1_grs, bins=20);\n",
    "axes[1].hist(y2_grs, bins=20);"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "ename": "AssertionError",
     "evalue": "Could not broadcast dimensions. Broadcasting is only allowed along axes that have a statically known length 1. Use `specify_broadcastable` to inform PyTensor of a known shape.\nApply node that caused the error: Assert{msg='Could not broadcast dimensi...'}(rating_reduced_dim, ScalarFromTensor.0)\nToposort index: 2\nInputs types: [TensorType(int64, shape=()), ScalarType(bool)]\nInputs shapes: [(), ()]\nInputs strides: [(), ()]\nInputs values: [array(3), np.False_]\nOutputs clients: [[MakeVector{dtype='int64'}(subject__factor_dim, treat__expr_dim, Assert{msg='Could not broadcast dimensi...'}.0)]]\n\nHINT: Re-running with most PyTensor optimizations disabled could provide a back-trace showing when this node was created. This can be done by setting the PyTensor flag 'optimizer=fast_compile'. If that does not work, PyTensor optimizations can be disabled with 'optimizer=None'.\nHINT: Use the PyTensor flag `exception_verbosity=high` for a debug print-out and storage map footprint of this Apply node.",
     "output_type": "error",
     "traceback": [
      "\u001b[31m---------------------------------------------------------------------------\u001b[39m",
      "\u001b[31mAssertionError\u001b[39m                            Traceback (most recent call last)",
      "\u001b[36mFile \u001b[39m\u001b[32m~/oss/bambinos/bambi/.pixi/envs/dev/lib/python3.13/site-packages/pytensor/compile/function/types.py:1039\u001b[39m, in \u001b[36mFunction.__call__\u001b[39m\u001b[34m(self, output_subset, *args, **kwargs)\u001b[39m\n\u001b[32m   1038\u001b[39m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[32m-> \u001b[39m\u001b[32m1039\u001b[39m     outputs = \u001b[43mvm\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m \u001b[38;5;28;01mif\u001b[39;00m output_subset \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m \u001b[38;5;28;01melse\u001b[39;00m vm(output_subset=output_subset)\n\u001b[32m   1040\u001b[39m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mException\u001b[39;00m:\n",
      "\u001b[31mAssertionError\u001b[39m: Could not broadcast dimensions. Broadcasting is only allowed along axes that have a statically known length 1. Use `specify_broadcastable` to inform PyTensor of a known shape.",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[31mAssertionError\u001b[39m                            Traceback (most recent call last)",
      "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[22]\u001b[39m\u001b[32m, line 8\u001b[39m\n\u001b[32m      3\u001b[39m inhaler[\u001b[33m\"\u001b[39m\u001b[33mtreat\u001b[39m\u001b[33m\"\u001b[39m] = pd.Categorical(inhaler[\u001b[33m\"\u001b[39m\u001b[33mtreat\u001b[39m\u001b[33m\"\u001b[39m])\n\u001b[32m      5\u001b[39m model = bmb.Model(\n\u001b[32m      6\u001b[39m     \u001b[33m\"\u001b[39m\u001b[33mrating ~ 1 + period + treat + (1 + treat|subject)\u001b[39m\u001b[33m\"\u001b[39m, inhaler, family=\u001b[33m\"\u001b[39m\u001b[33mcategorical\u001b[39m\u001b[33m\"\u001b[39m\n\u001b[32m      7\u001b[39m )\n\u001b[32m----> \u001b[39m\u001b[32m8\u001b[39m idata = \u001b[43mmodel\u001b[49m\u001b[43m.\u001b[49m\u001b[43mfit\u001b[49m\u001b[43m(\u001b[49m\u001b[43mtune\u001b[49m\u001b[43m=\u001b[49m\u001b[32;43m200\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdraws\u001b[49m\u001b[43m=\u001b[49m\u001b[32;43m200\u001b[39;49m\u001b[43m)\u001b[49m\n",
      "\u001b[36mFile \u001b[39m\u001b[32m~/oss/bambinos/bambi/bambi/models.py:327\u001b[39m, in \u001b[36mModel.fit\u001b[39m\u001b[34m(self, draws, tune, discard_tuned_samples, omit_offsets, include_mean, include_response_params, inference_method, init, n_init, chains, cores, random_seed, **kwargs)\u001b[39m\n\u001b[32m    324\u001b[39m         inference_method = method\n\u001b[32m    326\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28mself\u001b[39m.built:\n\u001b[32m--> \u001b[39m\u001b[32m327\u001b[39m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mbuild\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m    329\u001b[39m \u001b[38;5;66;03m# Tell user which event is being modeled\u001b[39;00m\n\u001b[32m    330\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(\u001b[38;5;28mself\u001b[39m.family, univariate.Bernoulli):\n",
      "\u001b[36mFile \u001b[39m\u001b[32m~/oss/bambinos/bambi/bambi/models.py:366\u001b[39m, in \u001b[36mModel.build\u001b[39m\u001b[34m(self)\u001b[39m\n\u001b[32m    361\u001b[39m \u001b[38;5;250m\u001b[39m\u001b[33;03m\"\"\"Set up the model for sampling/fitting\u001b[39;00m\n\u001b[32m    362\u001b[39m \n\u001b[32m    363\u001b[39m \u001b[33;03mCreates an instance of the underlying PyMC model and adds all the necessary terms to it.\u001b[39;00m\n\u001b[32m    364\u001b[39m \u001b[33;03m\"\"\"\u001b[39;00m\n\u001b[32m    365\u001b[39m \u001b[38;5;28mself\u001b[39m.backend = PyMCModel()\n\u001b[32m--> \u001b[39m\u001b[32m366\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mbackend\u001b[49m\u001b[43m.\u001b[49m\u001b[43mbuild\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[32m    367\u001b[39m \u001b[38;5;28mself\u001b[39m.built = \u001b[38;5;28;01mTrue\u001b[39;00m\n",
      "\u001b[36mFile \u001b[39m\u001b[32m~/oss/bambinos/bambi/bambi/backend/pymc.py:94\u001b[39m, in \u001b[36mPyMCModel.build\u001b[39m\u001b[34m(self, spec)\u001b[39m\n\u001b[32m     92\u001b[39m \u001b[38;5;28;01mfor\u001b[39;00m name, component \u001b[38;5;129;01min\u001b[39;00m spec.distributional_components.items():\n\u001b[32m     93\u001b[39m     \u001b[38;5;28mself\u001b[39m.components[name] = DistributionalComponent(component)\n\u001b[32m---> \u001b[39m\u001b[32m94\u001b[39m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mcomponents\u001b[49m\u001b[43m[\u001b[49m\u001b[43mname\u001b[49m\u001b[43m]\u001b[49m\u001b[43m.\u001b[49m\u001b[43mbuild\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mspec\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m     96\u001b[39m \u001b[38;5;66;03m# Add response\u001b[39;00m\n\u001b[32m     97\u001b[39m \u001b[38;5;28mself\u001b[39m.response_component = ResponseComponent(spec.response_component)\n",
      "\u001b[36mFile \u001b[39m\u001b[32m~/oss/bambinos/bambi/bambi/backend/model_components.py:57\u001b[39m, in \u001b[36mDistributionalComponent.build\u001b[39m\u001b[34m(self, pymc_backend, bmb_model)\u001b[39m\n\u001b[32m     55\u001b[39m \u001b[38;5;28mself\u001b[39m.build_common_terms(pymc_backend, bmb_model)\n\u001b[32m     56\u001b[39m \u001b[38;5;28mself\u001b[39m.build_hsgp_terms(bmb_model, pymc_backend)\n\u001b[32m---> \u001b[39m\u001b[32m57\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mbuild_group_specific_terms\u001b[49m\u001b[43m(\u001b[49m\u001b[43mpymc_backend\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mbmb_model\u001b[49m\u001b[43m)\u001b[49m\n",
      "\u001b[36mFile \u001b[39m\u001b[32m~/oss/bambinos/bambi/bambi/backend/model_components.py:141\u001b[39m, in \u001b[36mDistributionalComponent.build_group_specific_terms\u001b[39m\u001b[34m(self, pymc_backend, bmb_model)\u001b[39m\n\u001b[32m    138\u001b[39m         pymc_backend.model.add_coords({name: values})\n\u001b[32m    140\u001b[39m \u001b[38;5;66;03m# Build\u001b[39;00m\n\u001b[32m--> \u001b[39m\u001b[32m141\u001b[39m coef, predictor = \u001b[43mgroup_specific_term\u001b[49m\u001b[43m.\u001b[49m\u001b[43mbuild\u001b[49m\u001b[43m(\u001b[49m\u001b[43mbmb_model\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m    143\u001b[39m \u001b[38;5;66;03m# Add to the linear predictor\u001b[39;00m\n\u001b[32m    144\u001b[39m \u001b[38;5;66;03m# The loop through predictor columns is not the most beautiful alternative.\u001b[39;00m\n\u001b[32m    145\u001b[39m \u001b[38;5;66;03m# But it's the fastest. Doing matrix multiplication, pm.math.dot(data, coef), is slower.\u001b[39;00m\n\u001b[32m    146\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m predictor.ndim > \u001b[32m1\u001b[39m:\n",
      "\u001b[36mFile \u001b[39m\u001b[32m~/oss/bambinos/bambi/bambi/backend/terms.py:117\u001b[39m, in \u001b[36mGroupSpecificTerm.build\u001b[39m\u001b[34m(self, spec)\u001b[39m\n\u001b[32m    112\u001b[39m coef = \u001b[38;5;28mself\u001b[39m.build_distribution(\u001b[38;5;28mself\u001b[39m.term.prior, label, dims=dims, **kwargs)\n\u001b[32m    114\u001b[39m \u001b[38;5;66;03m# Squeeze ensures we don't have a shape of (n, 1) when we mean (n, )\u001b[39;00m\n\u001b[32m    115\u001b[39m \u001b[38;5;66;03m# This happens with categorical predictors with two levels and intercept.\u001b[39;00m\n\u001b[32m    116\u001b[39m \u001b[38;5;66;03m# See https://github.com/pymc-devs/pymc/issues/7246\u001b[39;00m\n\u001b[32m--> \u001b[39m\u001b[32m117\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(\u001b[43mcoef\u001b[49m\u001b[43m.\u001b[49m\u001b[43mshape\u001b[49m\u001b[43m.\u001b[49m\u001b[43meval\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m) == \u001b[32m2\u001b[39m \u001b[38;5;129;01mand\u001b[39;00m coef.shape.eval()[-\u001b[32m1\u001b[39m] == \u001b[32m1\u001b[39m:\n\u001b[32m    118\u001b[39m     coef = pt.specify_broadcastable(coef, \u001b[32m1\u001b[39m).squeeze()\n\u001b[32m    119\u001b[39m coef = coef[\u001b[38;5;28mself\u001b[39m.term.group_index]\n",
      "\u001b[36mFile \u001b[39m\u001b[32m~/oss/bambinos/bambi/.pixi/envs/dev/lib/python3.13/site-packages/pytensor/graph/basic.py:662\u001b[39m, in \u001b[36mVariable.eval\u001b[39m\u001b[34m(self, inputs_to_values, **kwargs)\u001b[39m\n\u001b[32m    656\u001b[39m         warnings.warn(\n\u001b[32m    657\u001b[39m             \u001b[33m\"\u001b[39m\u001b[33mKeyword arguments could not be used to create a cache key for the underlying variable. \u001b[39m\u001b[33m\"\u001b[39m\n\u001b[32m    658\u001b[39m             \u001b[33mf\u001b[39m\u001b[33m\"\u001b[39m\u001b[33mA function will be recompiled on every call with such keyword arguments.\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;132;01m{\u001b[39;00mexc\u001b[38;5;132;01m}\u001b[39;00m\u001b[33m\"\u001b[39m\n\u001b[32m    659\u001b[39m         )\n\u001b[32m    661\u001b[39m args = [parsed_inputs_to_values[param] \u001b[38;5;28;01mfor\u001b[39;00m param \u001b[38;5;129;01min\u001b[39;00m inputs]\n\u001b[32m--> \u001b[39m\u001b[32m662\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mfn\u001b[49m\u001b[43m(\u001b[49m\u001b[43m*\u001b[49m\u001b[43margs\u001b[49m\u001b[43m)\u001b[49m\n",
      "\u001b[36mFile \u001b[39m\u001b[32m~/oss/bambinos/bambi/.pixi/envs/dev/lib/python3.13/site-packages/pytensor/compile/function/types.py:1049\u001b[39m, in \u001b[36mFunction.__call__\u001b[39m\u001b[34m(self, output_subset, *args, **kwargs)\u001b[39m\n\u001b[32m   1047\u001b[39m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mhasattr\u001b[39m(\u001b[38;5;28mself\u001b[39m.vm, \u001b[33m\"\u001b[39m\u001b[33mthunks\u001b[39m\u001b[33m\"\u001b[39m):\n\u001b[32m   1048\u001b[39m         thunk = \u001b[38;5;28mself\u001b[39m.vm.thunks[\u001b[38;5;28mself\u001b[39m.vm.position_of_error]\n\u001b[32m-> \u001b[39m\u001b[32m1049\u001b[39m     \u001b[43mraise_with_op\u001b[49m\u001b[43m(\u001b[49m\n\u001b[32m   1050\u001b[39m \u001b[43m        \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mmaker\u001b[49m\u001b[43m.\u001b[49m\u001b[43mfgraph\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m   1051\u001b[39m \u001b[43m        \u001b[49m\u001b[43mnode\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mvm\u001b[49m\u001b[43m.\u001b[49m\u001b[43mnodes\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mvm\u001b[49m\u001b[43m.\u001b[49m\u001b[43mposition_of_error\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m   1052\u001b[39m \u001b[43m        \u001b[49m\u001b[43mthunk\u001b[49m\u001b[43m=\u001b[49m\u001b[43mthunk\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m   1053\u001b[39m \u001b[43m        \u001b[49m\u001b[43mstorage_map\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43mgetattr\u001b[39;49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mvm\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[33;43m\"\u001b[39;49m\u001b[33;43mstorage_map\u001b[39;49m\u001b[33;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mNone\u001b[39;49;00m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m   1054\u001b[39m \u001b[43m    \u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m   1055\u001b[39m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[32m   1056\u001b[39m     \u001b[38;5;66;03m# old-style linkers raise their own exceptions\u001b[39;00m\n\u001b[32m   1057\u001b[39m     \u001b[38;5;28;01mraise\u001b[39;00m\n",
      "\u001b[36mFile \u001b[39m\u001b[32m~/oss/bambinos/bambi/.pixi/envs/dev/lib/python3.13/site-packages/pytensor/link/utils.py:526\u001b[39m, in \u001b[36mraise_with_op\u001b[39m\u001b[34m(fgraph, node, thunk, exc_info, storage_map)\u001b[39m\n\u001b[32m    521\u001b[39m     warnings.warn(\n\u001b[32m    522\u001b[39m         \u001b[33mf\u001b[39m\u001b[33m\"\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mexc_type\u001b[38;5;132;01m}\u001b[39;00m\u001b[33m error does not allow us to add an extra error message\u001b[39m\u001b[33m\"\u001b[39m\n\u001b[32m    523\u001b[39m     )\n\u001b[32m    524\u001b[39m     \u001b[38;5;66;03m# Some exception need extra parameter in inputs. So forget the\u001b[39;00m\n\u001b[32m    525\u001b[39m     \u001b[38;5;66;03m# extra long error message in that case.\u001b[39;00m\n\u001b[32m--> \u001b[39m\u001b[32m526\u001b[39m \u001b[38;5;28;01mraise\u001b[39;00m exc_value.with_traceback(exc_trace)\n",
      "\u001b[36mFile \u001b[39m\u001b[32m~/oss/bambinos/bambi/.pixi/envs/dev/lib/python3.13/site-packages/pytensor/compile/function/types.py:1039\u001b[39m, in \u001b[36mFunction.__call__\u001b[39m\u001b[34m(self, output_subset, *args, **kwargs)\u001b[39m\n\u001b[32m   1037\u001b[39m     t0_fn = time.perf_counter()\n\u001b[32m   1038\u001b[39m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[32m-> \u001b[39m\u001b[32m1039\u001b[39m     outputs = \u001b[43mvm\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m \u001b[38;5;28;01mif\u001b[39;00m output_subset \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m \u001b[38;5;28;01melse\u001b[39;00m vm(output_subset=output_subset)\n\u001b[32m   1040\u001b[39m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mException\u001b[39;00m:\n\u001b[32m   1041\u001b[39m     \u001b[38;5;28mself\u001b[39m._restore_defaults()\n",
      "\u001b[31mAssertionError\u001b[39m: Could not broadcast dimensions. Broadcasting is only allowed along axes that have a statically known length 1. Use `specify_broadcastable` to inform PyTensor of a known shape.\nApply node that caused the error: Assert{msg='Could not broadcast dimensi...'}(rating_reduced_dim, ScalarFromTensor.0)\nToposort index: 2\nInputs types: [TensorType(int64, shape=()), ScalarType(bool)]\nInputs shapes: [(), ()]\nInputs strides: [(), ()]\nInputs values: [array(3), np.False_]\nOutputs clients: [[MakeVector{dtype='int64'}(subject__factor_dim, treat__expr_dim, Assert{msg='Could not broadcast dimensi...'}.0)]]\n\nHINT: Re-running with most PyTensor optimizations disabled could provide a back-trace showing when this node was created. This can be done by setting the PyTensor flag 'optimizer=fast_compile'. If that does not work, PyTensor optimizations can be disabled with 'optimizer=None'.\nHINT: Use the PyTensor flag `exception_verbosity=high` for a debug print-out and storage map footprint of this Apply node."
     ]
    }
   ],
   "source": [
    "inhaler = pd.read_csv(\"../../tests/data/inhaler.csv\")\n",
    "inhaler[\"rating\"] = pd.Categorical(inhaler[\"rating\"], categories=[1, 2, 3, 4])\n",
    "inhaler[\"treat\"] = pd.Categorical(inhaler[\"treat\"])\n",
    "\n",
    "model = bmb.Model(\n",
    "    \"rating ~ 1 + period + treat + (1 + treat|subject)\", inhaler, family=\"categorical\"\n",
    ")\n",
    "idata = model.fit(tune=200, draws=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>subject</th>\n",
       "      <th>rating</th>\n",
       "      <th>treat</th>\n",
       "      <th>period</th>\n",
       "      <th>carry</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.5</td>\n",
       "      <td>0.5</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>999</td>\n",
       "      <td>1</td>\n",
       "      <td>0.5</td>\n",
       "      <td>0.5</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   subject rating treat  period  carry\n",
       "0        1      1   0.5     0.5      0\n",
       "1      999      1   0.5     0.5      0"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_new = inhaler.head(2).reset_index(drop=True)\n",
    "df_new[\"subject\"] = [1, 999]\n",
    "df_new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = model.predict(idata, data=df_new, inplace=False, sample_new_groups=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA9oAAALgCAYAAABrgQ84AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABxvklEQVR4nO3de3xU1b3//3cuZJJALgRMJtGAXJT7RUHj1BuXlBD4Uj3knIoixsoBtcFW0ipNiyighqJV1BPhtAcBv0Kxnh/QFjkgoECVgBjNEcGmBtGgMEHEZEiAyW3//vCbqSPXSdbM5PJ6Ph77UWbvNWvWXo355D37FmJZliUAAAAAAGBEaLAHAAAAAABAW0LQBgAAAADAIII2AAAAAAAGEbQBAAAAADCIoA0AAAAAgEEEbQAAAAAADCJoAwAAAABgUHiwB9AUDQ0NOnz4sGJiYhQSEhLs4QAAIMuydOLECaWkpCg0lO+xm4taDwBoaXyp9a0yaB8+fFipqanBHgYAAGc4dOiQLrvssmAPo9Wj1gMAWqqLqfWtMmjHxMRI+nYHY2NjgzwaAAAkl8ul1NRUT41C81DrAQAtjS+1vlUG7cZTyGJjYym+AIAWhdOczaDWAwBaqoup9VxEBgAAAACAQQRtAAAAAAAMImgDAAAAAGAQQRsAAAAAAIMI2gAAAAAAGETQBgAAAADAIII2AAAAAAAGEbQBAAAAADCIoA0AAAAAgEEEbQAAAAAADCJoAwAAAABgUHiwBwAAwMX6suKUvqmuMdpn544RujQ+ymifaJlM//zwswMAOBeCNgCgVfiy4pTSf7ddp2rrjfYb1SFMW35xM4GpjfPHzw8/OwCAcyFoAwBahW+qa3Sqtl6Lbhuq3omdjPRZerRKD75arG+qawhLbZzpnx9+dgAA50PQBgC0Kr0TO2ngpXHBHgZaKX5+AACB4PPN0Hbs2KEJEyYoJSVFISEhWrdundf2kJCQsy5PPfWUp83ll19+xvYFCxY0e2cAAAAAAAg2n4N2dXW1hgwZooKCgrNuP3LkiNfy0ksvKSQkRFlZWV7t5s2b59XugQceaNoeAAAAAADQgvh86nhmZqYyMzPPud1ut3u9/vOf/6yRI0eqZ8+eXutjYmLOaAsAAAAAQGvn1+dol5eX6/XXX9fUqVPP2LZgwQJ16dJFV111lZ566inV1dWdsx+32y2Xy+W1AAAAAADQEvn1ZmgrVqxQTEyMJk6c6LX+Zz/7ma6++molJCRo586dysvL05EjR/TMM8+ctZ/8/HzNnTvXn0MFAAAAAMAIvwbtl156SZMnT1ZkZKTX+tzcXM+/Bw8erIiICN17773Kz8+XzWY7o5+8vDyv97hcLqWmpvpv4AAAAAAANJHfgvbf/vY3lZSU6NVXX71g27S0NNXV1emzzz5Tnz59zthus9nOGsABAAAAAGhp/HaN9tKlSzVs2DANGTLkgm2Li4sVGhqqxMREfw0HAAAAAICA8PmIdlVVlUpLSz2vDx48qOLiYiUkJKhbt26Svj21+7XXXtPvfve7M95fWFio3bt3a+TIkYqJiVFhYaFmzpypO++8U507d27GrgAAAAAAEHw+B+333ntPI0eO9LxuvHY6Oztby5cvlyStXr1almXp9ttvP+P9NptNq1ev1mOPPSa3260ePXpo5syZXtdgAwAAAADQWvkctEeMGCHLss7bZvr06Zo+ffpZt1199dXatWuXrx8LAAAAAECr4NfnaAMAAAAA0N4QtAEAAAAAMIigDQAAPPLz83XNNdcoJiZGiYmJuvXWW1VSUuLV5vTp08rJyVGXLl3UqVMnZWVlqby83KtNWVmZxo8fr+joaCUmJuqhhx5SXV1dIHcFAICgIWgDAACP7du3KycnR7t27dLmzZtVW1urMWPGqLq62tNm5syZ+utf/6rXXntN27dv1+HDhzVx4kTP9vr6eo0fP141NTXauXOnVqxYoeXLl2vOnDnB2CUAAALO55uhAQCAtmvjxo1er5cvX67ExEQVFRXppptuUmVlpZYuXapVq1Zp1KhRkqRly5apX79+2rVrl6677jq98cYb2r9/v7Zs2aKkpCQNHTpU8+fP16xZs/TYY48pIiIiGLsGAEDAcEQbAACcU2VlpSQpISFBklRUVKTa2lqlp6d72vTt21fdunVTYWGhJKmwsFCDBg1SUlKSp01GRoZcLpf27dt31s9xu91yuVxeCwAArRVBGwAAnFVDQ4MefPBBXX/99Ro4cKAkyel0KiIiQvHx8V5tk5KS5HQ6PW2+G7IbtzduO5v8/HzFxcV5ltTUVMN7AwBA4BC0AQDAWeXk5Oijjz7S6tWr/f5ZeXl5qqys9CyHDh3y+2cCAOAvXKMNAADOMGPGDK1fv147duzQZZdd5llvt9tVU1OjiooKr6Pa5eXlstvtnjbvvvuuV3+NdyVvbPN9NptNNpvN8F4AABAcHNEGAAAelmVpxowZWrt2rd5880316NHDa/uwYcPUoUMHbd261bOupKREZWVlcjgckiSHw6G9e/fq6NGjnjabN29WbGys+vfvH5gdAQAgiDiiDQAAPHJycrRq1Sr9+c9/VkxMjOea6ri4OEVFRSkuLk5Tp05Vbm6uEhISFBsbqwceeEAOh0PXXXedJGnMmDHq37+/pkyZooULF8rpdGr27NnKycnhqDUAoF0gaAMAAI/FixdLkkaMGOG1ftmyZbr77rslSc8++6xCQ0OVlZUlt9utjIwMvfjii562YWFhWr9+ve6//345HA517NhR2dnZmjdvXqB2AwCAoCJoAwAAD8uyLtgmMjJSBQUFKigoOGeb7t27a8OGDSaHBgBAq8E12gAAAAAAGETQBgAAAADAIII2AAAAAAAGEbQBAAAAADCIoA0AAAAAgEEEbQAAAAAADCJoAwAAAABgEEEbAAAAAACDCNoAAAAAABhE0AYAAAAAwCCCNgAAAAAABhG0AQAAAAAwiKANAAAAAIBBBG0AAAAAAAwiaAMAAAAAYBBBGwAAAAAAgwjaAAAAAAAYRNAGAAAAAMAgn4P2jh07NGHCBKWkpCgkJETr1q3z2n733XcrJCTEaxk7dqxXm+PHj2vy5MmKjY1VfHy8pk6dqqqqqmbtCAAAAAAALYHPQbu6ulpDhgxRQUHBOduMHTtWR44c8Sx//OMfvbZPnjxZ+/bt0+bNm7V+/Xrt2LFD06dP9330AAAAAAC0MOG+viEzM1OZmZnnbWOz2WS328+67eOPP9bGjRu1Z88eDR8+XJL0wgsvaNy4cXr66aeVkpLi65AAAAAAAGgx/HKN9rZt25SYmKg+ffro/vvv19dff+3ZVlhYqPj4eE/IlqT09HSFhoZq9+7dZ+3P7XbL5XJ5LQAAAAAAtETGg/bYsWP18ssva+vWrfrtb3+r7du3KzMzU/X19ZIkp9OpxMREr/eEh4crISFBTqfzrH3m5+crLi7Os6SmppoeNgAAAAAARvh86viFTJo0yfPvQYMGafDgwerVq5e2bdum0aNHN6nPvLw85ebmel67XC7CNgAAAACgRfL747169uyprl27qrS0VJJkt9t19OhRrzZ1dXU6fvz4Oa/rttlsio2N9VoAAIB/XOgJI99/ukjj8tRTT3naXH755WdsX7BgQYD3BACA4PB70P7iiy/09ddfKzk5WZLkcDhUUVGhoqIiT5s333xTDQ0NSktL8/dwAADABVzoCSPffbLIkSNH9NJLLykkJERZWVle7ebNm+fV7oEHHgjE8AEACDqfTx2vqqryHJ2WpIMHD6q4uFgJCQlKSEjQ3LlzlZWVJbvdrgMHDujhhx9W7969lZGRIUnq16+fxo4dq2nTpmnJkiWqra3VjBkzNGnSJO44DgBAC3ChJ4x8/wy0P//5zxo5cqR69uzptT4mJuacZ6sBANCW+XxE+7333tNVV12lq666SpKUm5urq666SnPmzFFYWJg+/PBD/ehHP9KVV16pqVOnatiwYfrb3/4mm83m6WPlypXq27evRo8erXHjxumGG27Q73//e3N7BQAAAqK8vFyvv/66pk6desa2BQsWqEuXLrrqqqv01FNPqa6u7pz98IQRAEBb4vMR7REjRsiyrHNu37Rp0wX7SEhI0KpVq3z9aAAA0MKsWLFCMTExmjhxotf6n/3sZ7r66quVkJCgnTt3Ki8vT0eOHNEzzzxz1n7y8/M1d+7cQAwZAAC/M37XcQAA0H689NJLmjx5siIjI73Wf/dpIYMHD1ZERITuvfde5efne53l1ognjAAA2hKCNgAAaJK//e1vKikp0auvvnrBtmlpaaqrq9Nnn32mPn36nLHdZrOdNYADANAa+f2u4wAAoG1aunSphg0bpiFDhlywbXFxsUJDQ5WYmBiAkQEAEFwc0QYAAF7O94SRbt26Sfr21O7XXntNv/vd7854f2FhoXbv3q2RI0cqJiZGhYWFmjlzpu6880517tw5YPsBAECwELQBAICX9957TyNHjvS8brx2Ojs7W8uXL5ckrV69WpZl6fbbbz/j/TabTatXr9Zjjz0mt9utHj16aObMmV7XYAMA0JYRtAEAgJcLPWFEkqZPn67p06efddvVV1+tXbt2+WNoAAC0ClyjDQAAAACAQQRtAAAAAAAMImgDAAAAAGAQQRsAAAAAAIMI2gAAAAAAGETQBgAAAADAIII2AAAAAAAGEbQBAAAAADCIoA0AAAAAgEEEbQAAAAAADCJoAwAAAABgEEEbAAAAAACDCNoAAAAAABhE0AYAAAAAwCCCNgAAAAAABhG0AQAAAAAwiKANAAAAAIBBBG0AAAAAAAwiaAMAAAAAYBBBGwAAAAAAgwjaAAAAAAAYRNAGAAAAAMAggjYAAAAAAAYRtAEAAAAAMIigDQAAvOzYsUMTJkxQSkqKQkJCtG7dOq/td999t0JCQryWsWPHerU5fvy4Jk+erNjYWMXHx2vq1KmqqqoK4F4AABA8Pgft8xXf2tpazZo1S4MGDVLHjh2VkpKiu+66S4cPH/bq4/LLLz+jQC9YsKDZOwMAAJqvurpaQ4YMUUFBwTnbjB07VkeOHPEsf/zjH722T548Wfv27dPmzZu1fv167dixQ9OnT/f30AEAaBHCfX1DY/G95557NHHiRK9tJ0+e1Pvvv69HHnlEQ4YM0TfffKOf//zn+tGPfqT33nvPq+28efM0bdo0z+uYmJgm7gIAADApMzNTmZmZ521js9lkt9vPuu3jjz/Wxo0btWfPHg0fPlyS9MILL2jcuHF6+umnlZKSYnzMAAC0JD4H7fMV37i4OG3evNlr3X/8x3/o2muvVVlZmbp16+ZZHxMTc84CDQAAWrZt27YpMTFRnTt31qhRo/T444+rS5cukqTCwkLFx8d7QrYkpaenKzQ0VLt379a//Mu/nNGf2+2W2+32vHa5XP7fCQAA/MTv12hXVlYqJCRE8fHxXusXLFigLl266KqrrtJTTz2lurq6c/bhdrvlcrm8FgAAEBxjx47Vyy+/rK1bt+q3v/2ttm/frszMTNXX10uSnE6nEhMTvd4THh6uhIQEOZ3Os/aZn5+vuLg4z5Kamur3/QAAwF98PqLti9OnT2vWrFm6/fbbFRsb61n/s5/9TFdffbUSEhK0c+dO5eXl6ciRI3rmmWfO2k9+fr7mzp3rz6ECAICLNGnSJM+/Bw0apMGDB6tXr17atm2bRo8e3aQ+8/LylJub63ntcrkI2wCAVstvQbu2tlY//vGPZVmWFi9e7LXtu4V08ODBioiI0L333qv8/HzZbLYz+qL4AgDQcvXs2VNdu3ZVaWmpRo8eLbvdrqNHj3q1qaur0/Hjx8952ZjNZjvr3wAAALRGfjl1vDFkf/7559q8ebPX0eyzSUtLU11dnT777LOzbrfZbIqNjfVaAABAy/DFF1/o66+/VnJysiTJ4XCooqJCRUVFnjZvvvmmGhoalJaWFqxhAgAQMMaPaDeG7E8++URvvfWW58Yo51NcXKzQ0NAzrucCAACBV1VVpdLSUs/rgwcPqri4WAkJCUpISNDcuXOVlZUlu92uAwcO6OGHH1bv3r2VkZEhSerXr5/Gjh2radOmacmSJaqtrdWMGTM0adIk7jgOAGgXfA7a5yu+ycnJ+td//Ve9//77Wr9+verr6z03PUlISFBERIQKCwu1e/dujRw5UjExMSosLNTMmTN15513qnPnzub2DAAANMl7772nkSNHel43Xr6VnZ2txYsX68MPP9SKFStUUVGhlJQUjRkzRvPnz/c69XvlypWaMWOGRo8erdDQUGVlZen5558P+L4AABAMPgft8xXfxx57TH/5y18kSUOHDvV631tvvaURI0bIZrNp9erVeuyxx+R2u9WjRw/NnDnT6xpsAAAQPCNGjJBlWefcvmnTpgv2kZCQoFWrVpkcFgAArYbPQftCxfd82yTp6quv1q5du3z9WAAAAAAAWgW/P0cbAAAAAID2hKANAAAAAIBBBG0AAAAAAAwiaAMAAAAAYBBBGwAAAAAAgwjaAAAAAAAYRNAGAAAAAMAggjYAAAAAAAYRtAEAAAAAMIigDQAAAACAQQRtAAAAAAAMImgDAAAAAGAQQRsAAAAAAIMI2gAAAAAAGETQBgAAAADAIII2AAAAAAAGEbQBAAAAADCIoA0AAAAAgEEEbQAAAAAADCJoAwAAAABgEEEbAAAAAACDCNoAAAAAABhE0AYAAAAAwCCCNgAA8LJjxw5NmDBBKSkpCgkJ0bp16zzbamtrNWvWLA0aNEgdO3ZUSkqK7rrrLh0+fNirj8svv1whISFey4IFCwK8JwAABAdBGwAAeKmurtaQIUNUUFBwxraTJ0/q/fff1yOPPKL3339fa9asUUlJiX70ox+d0XbevHk6cuSIZ3nggQcCMXwAAIIuPNgDAAAALUtmZqYyMzPPui0uLk6bN2/2Wvcf//Efuvbaa1VWVqZu3bp51sfExMhut/t1rAAAtEQc0QYAAM1SWVmpkJAQxcfHe61fsGCBunTpoquuukpPPfWU6urqztmH2+2Wy+XyWgAAaK04og0AAJrs9OnTmjVrlm6//XbFxsZ61v/sZz/T1VdfrYSEBO3cuVN5eXk6cuSInnnmmbP2k5+fr7lz5wZq2AAA+BVBGwAANEltba1+/OMfy7IsLV682Gtbbm6u59+DBw9WRESE7r33XuXn58tms53RV15entd7XC6XUlNT/Td4AAD8iKANAAB81hiyP//8c7355pteR7PPJi0tTXV1dfrss8/Up0+fM7bbbLazBnAAAFojgjYAAPBJY8j+5JNP9NZbb6lLly4XfE9xcbFCQ0OVmJgYgBECABBcPt8M7XzP1pQky7I0Z84cJScnKyoqSunp6frkk0+82hw/flyTJ09WbGys4uPjNXXqVFVVVTVrRwAAgBlVVVUqLi5WcXGxJOngwYMqLi5WWVmZamtr9a//+q967733tHLlStXX18vpdMrpdKqmpkaSVFhYqEWLFul///d/9emnn2rlypWaOXOm7rzzTnXu3DmIewYAQGD4HLTP92xNSVq4cKGef/55LVmyRLt371bHjh2VkZGh06dPe9pMnjxZ+/bt0+bNm7V+/Xrt2LFD06dPb/peAAAAY9577z1dddVVuuqqqyR9e731VVddpTlz5ujLL7/UX/7yF33xxRcaOnSokpOTPcvOnTslfXsa+OrVq3XzzTdrwIABeuKJJzRz5kz9/ve/D+ZuAQAQMD6fOn6+Z2talqVFixZp9uzZuuWWWyRJL7/8spKSkrRu3TpNmjRJH3/8sTZu3Kg9e/Zo+PDhkqQXXnhB48aN09NPP62UlJRm7A4AAGiuESNGyLKsc24/3zZJuvrqq7Vr1y7Tw2q2FB1T5LG9UkinZvcVeaxKKTpmYFQAgLbI6DXaBw8elNPpVHp6umddXFyc0tLSVFhYqEmTJqmwsFDx8fGekC1J6enpCg0N1e7du/Uv//IvZ/Trdrvldrs9r3m2JgAA8EWHqi+1xfaQote6L9z4IvSWtMVm06GqayTFGekTANB2GA3aTqdTkpSUlOS1PikpybPN6XSecSOU8PBwJSQkeNp8H8/WBAAAzRF2+riiQ9w6NPI5pV4xtNn9HfqkWKlv/Vxhp483f3AAgDanVdx1nGdrAgAAE9zxvaWUoc3v5ytu4goAODefb4Z2Pna7XZJUXl7utb68vNyzzW636+jRo17b6+rqdPz4cU+b77PZbIqNjfVaAAAAAABoiYwG7R49eshut2vr1q2edS6XS7t375bD4ZAkORwOVVRUqKioyNPmzTffVENDg9LS0kwOBwAAAACAgPP51PGqqiqVlpZ6Xjc+WzMhIUHdunXTgw8+qMcff1xXXHGFevTooUceeUQpKSm69dZbJUn9+vXT2LFjNW3aNC1ZskS1tbWaMWOGJk2axB3HAQAAAACtns9B+7333tPIkSM9rxuvnc7Oztby5cv18MMPq7q6WtOnT1dFRYVuuOEGbdy4UZGRkZ73rFy5UjNmzNDo0aMVGhqqrKwsPf/88wZ2BwAAAACA4PI5aF/o2ZohISGaN2+e5s2bd842CQkJWrVqla8fDQAAAABAi2f0Gm0AAAAAANo7gjYAAAAAAAYRtAEAAAAAMIigDQAAAACAQQRtAAAAAAAMImgDAAAAAGAQQRsAAAAAAIMI2gAAAAAAGBQe7AEAAHCxUnRMkcf2SiGdjPQXeaxKKTpmpC8AAIBGBG0AQKvQoepLbbE9pOi1bmN99pa0xWbToaprJMUZ6xcAALRvBG0AQKsQdvq4okPcOjTyOaVeMdRIn4c+KVbqWz9X2OnjRvoDAACQCNoAgFbGHd9bShlqpq+vqoz0AwAA8F3cDA0AAAAAAIMI2gAAAAAAGETQBgAAAADAIII2AAAAAAAGEbQBAICXHTt2aMKECUpJSVFISIjWrVvntd2yLM2ZM0fJycmKiopSenq6PvnkE682x48f1+TJkxUbG6v4+HhNnTpVVVXcfA4A0D4QtAEAgJfq6moNGTJEBQUFZ92+cOFCPf/881qyZIl2796tjh07KiMjQ6dPn/a0mTx5svbt26fNmzdr/fr12rFjh6ZPnx6oXQAAIKh4vBcAAPCSmZmpzMzMs26zLEuLFi3S7Nmzdcstt0iSXn75ZSUlJWndunWaNGmSPv74Y23cuFF79uzR8OHDJUkvvPCCxo0bp6efflopKSkB2xcAAIKBI9oAAOCiHTx4UE6nU+np6Z51cXFxSktLU2FhoSSpsLBQ8fHxnpAtSenp6QoNDdXu3bvP2q/b7ZbL5fJaAABorQjaAADgojmdTklSUlKS1/qkpCTPNqfTqcTERK/t4eHhSkhI8LT5vvz8fMXFxXmW1NRUP4weAIDAIGgDAICgy8vLU2VlpWc5dOhQsIcEAECTEbQBAMBFs9vtkqTy8nKv9eXl5Z5tdrtdR48e9dpeV1en48ePe9p8n81mU2xsrNcCAEBrRdAGAAAXrUePHrLb7dq6datnncvl0u7du+VwOCRJDodDFRUVKioq8rR588031dDQoLS0tICPGQCAQOOu4wAAwEtVVZVKS0s9rw8ePKji4mIlJCSoW7duevDBB/X444/riiuuUI8ePfTII48oJSVFt956qySpX79+Gjt2rKZNm6YlS5aotrZWM2bM0KRJk7jjOACgXSBoAwAAL++9955GjhzpeZ2bmytJys7O1vLly/Xwww+rurpa06dPV0VFhW644QZt3LhRkZGRnvesXLlSM2bM0OjRoxUaGqqsrCw9//zzAd8XAACCgaANAAC8jBgxQpZlnXN7SEiI5s2bp3nz5p2zTUJCglatWuWP4QEA0OJxjTYAAAAAAAYRtAEAAAAAMIigDQAAAACAQQRtAAAAAAAMMh60L7/8coWEhJyx5OTkSPr2Bivf33bfffeZHgYAAAAAAEFh/K7je/bsUX19vef1Rx99pB/+8If6t3/7N8+6adOmed2pNDo62vQwAAAAAAAICuNB+5JLLvF6vWDBAvXq1Us333yzZ110dLTsdrvpjwYAAAAAIOj8eo12TU2NXnnlFd1zzz0KCQnxrF+5cqW6du2qgQMHKi8vTydPnjxvP263Wy6Xy2sBAAAAAKAlMn5E+7vWrVuniooK3X333Z51d9xxh7p3766UlBR9+OGHmjVrlkpKSrRmzZpz9pOfn6+5c+f6c6gAAAAAABjh16C9dOlSZWZmKiUlxbNu+vTpnn8PGjRIycnJGj16tA4cOKBevXqdtZ+8vDzl5uZ6XrtcLqWmpvpv4AAAAAAANJHfgvbnn3+uLVu2nPdItSSlpaVJkkpLS88ZtG02m2w2m/ExAgAAAABgmt+u0V62bJkSExM1fvz487YrLi6WJCUnJ/trKAAAAAAABIxfjmg3NDRo2bJlys7OVnj4Pz/iwIEDWrVqlcaNG6cuXbroww8/1MyZM3XTTTdp8ODB/hgKAAAAAAAB5ZegvWXLFpWVlemee+7xWh8REaEtW7Zo0aJFqq6uVmpqqrKysjR79mx/DAMAAAAAgIDzS9AeM2aMLMs6Y31qaqq2b9/uj48EAAAAAKBF8OtztAEAAAAAaG8I2gAAAAAAGETQBgAAAADAIII2AAAAAAAGEbQBAAAAADDIL3cdBwBAkr6sOKVvqmuM9PXV8VPqbaQnAAAA/yJoAwD84suKU0r/3Xadqq030t+AkIMaaZNiozoY6Q8AAMBfCNoAAL/4prpGp2rrtei2oeqd2KnZ/UUei5PWSomdbAZGBwAA4D8EbQCAX/VO7KSBl8Y1v6OQ5od1AACAQOBmaAAAwCeXX365QkJCzlhycnIkSSNGjDhj23333RfkUQMAEDgc0QYAAD7Zs2eP6uv/ee39Rx99pB/+8If6t3/7N8+6adOmad68eZ7X0dHRAR0jAADBRNAGAAA+ueSSS7xeL1iwQL169dLNN9/sWRcdHS273R7ooQEA0CJw6jgAAGiympoavfLKK7rnnnsUEhLiWb9y5Up17dpVAwcOVF5enk6ePHneftxut1wul9cCAEBrxRFtAADQZOvWrVNFRYXuvvtuz7o77rhD3bt3V0pKij788EPNmjVLJSUlWrNmzTn7yc/P19y5cwMwYgAA/I+gDQAAmmzp0qXKzMxUSkqKZ9306dM9/x40aJCSk5M1evRoHThwQL169TprP3l5ecrNzfW8drlcSk1N9d/AAQDwI4I2AABoks8//1xbtmw575FqSUpLS5MklZaWnjNo22w22Ww8Ix0A0DZwjTYAAGiSZcuWKTExUePHjz9vu+LiYklScnJyAEYFAEDwcUQbAAD4rKGhQcuWLVN2drbCw//558SBAwe0atUqjRs3Tl26dNGHH36omTNn6qabbtLgwYODOGIAAAKHoA0AAHy2ZcsWlZWV6Z577vFaHxERoS1btmjRokWqrq5WamqqsrKyNHv27CCNFACAwCNoAwAAn40ZM0aWZZ2xPjU1Vdu3bw/CiAAAaDm4RhsAAAAAAIMI2gAAAAAAGETQBgAAAADAIII2AAAAAAAGEbQBAAAAADCIoA0AAAAAgEEEbQAAAAAADCJoAwAAAABgEEEbAAAAAACDCNoAAAAAABgUHuwBAAAAtFaHjp/S6S8rjfTVuWOELo2PMtIXACC4jAftxx57THPnzvVa16dPH/3973+XJJ0+fVq/+MUvtHr1arndbmVkZOjFF19UUlKS6aEAAAD4RWxUB0nS02+UaN+mGiN9RnUI05Zf3EzYBoA2wC9HtAcMGKAtW7b880PC//kxM2fO1Ouvv67XXntNcXFxmjFjhiZOnKh33nnHH0MBAAAwLrGTTZL03KShOt11ULP7Kz1apQdfLdY31TUEbQBoA/wStMPDw2W3289YX1lZqaVLl2rVqlUaNWqUJGnZsmXq16+fdu3apeuuu+6s/bndbrndbs9rl8vlj2EDAAD4pPclnaSUuGAPAwDQwvjlZmiffPKJUlJS1LNnT02ePFllZWWSpKKiItXW1io9Pd3Ttm/fvurWrZsKCwvP2V9+fr7i4uI8S2pqqj+GDQAAAABAsxkP2mlpaVq+fLk2btyoxYsX6+DBg7rxxht14sQJOZ1ORUREKD4+3us9SUlJcjqd5+wzLy9PlZWVnuXQoUOmhw0AAAAAgBHGTx3PzMz0/Hvw4MFKS0tT9+7d9ac//UlRUU275shms8lms5kaIgAAAAAAfuP352jHx8fryiuvVGlpqex2u2pqalRRUeHVpry8/KzXdAMAAAAA0Nr4PWhXVVXpwIEDSk5O1rBhw9ShQwdt3brVs72kpERlZWVyOBz+HgoAAAAAAH5n/NTxX/7yl5owYYK6d++uw4cP69FHH1VYWJhuv/12xcXFaerUqcrNzVVCQoJiY2P1wAMPyOFwnPOO4wCA1itFxxR5bK8U0qn5nR37R/P7AAAACADjQfuLL77Q7bffrq+//lqXXHKJbrjhBu3atUuXXHKJJOnZZ59VaGiosrKy5Ha7lZGRoRdffNH0MAAAQdah6kttsT2k6LXuCze+6E6jpegu5voDAADwA+NBe/Xq1efdHhkZqYKCAhUUFJj+aABACxJ2+riiQ9w6NPI5pV4x1Eyn0V2keB7xCAAAWjbjQRsAgO9yx/eWUoYGexgAAAAB4/eboQEAAAAA0J4QtAEAgE8ee+wxhYSEeC19+/b1bD99+rRycnLUpUsXderUSVlZWSovLw/iiAEACCxOHQcA/FPFIenk10a6slWUGukHLdOAAQO0ZcsWz+vw8H/+STFz5ky9/vrreu211xQXF6cZM2Zo4sSJeuedd4IxVAAAAo6gDQD4VsUhqeBaqfakke5SJZ20bKqPTDDSH1qW8PBw2e32M9ZXVlZq6dKlWrVqlUaNGiVJWrZsmfr166ddu3bxOE8AQLtA0AYAfOvk19+G7Il/kLpe2ezuSr+q0l1/PKDfd7rUwODQ0nzyySdKSUlRZGSkHA6H8vPz1a1bNxUVFam2tlbp6emetn379lW3bt1UWFh4zqDtdrvldv/zUXAul8vv+wAAgL8QtAEA3rpeaeQu4aetSh1WZfPHgxYnLS1Ny5cvV58+fXTkyBHNnTtXN954oz766CM5nU5FREQoPj7e6z1JSUlyOp3n7DM/P19z587188gBAAgMgjYAAPBJZmam59+DBw9WWlqaunfvrj/96U+KiopqUp95eXnKzc31vHa5XEpN5ZnpAIDWiaANAK3YlxWn9E11jZG+Io9VqbeRntDexMfH68orr1Rpaal++MMfqqamRhUVFV5HtcvLy896TXcjm80mm80WgNECAOB/BG0AaKW+rDil9N9t16naeiP9DQg5qNdt0tEqtxKN9Ij2oqqqSgcOHNCUKVM0bNgwdejQQVu3blVWVpYkqaSkRGVlZXI4HEEeKQAAgUHQBoBW6pvqGp2qrdei24aqd2KnZvf31T8ipO3Svi9dOtqx+ddWlx6tanYfaJl++ctfasKECerevbsOHz6sRx99VGFhYbr99tsVFxenqVOnKjc3VwkJCYqNjdUDDzwgh8PBHccBAO0GQRsAWrneiZ008NK4ZvdztDpWkvT0GyXat8nM6ehRHcLUuWOEkb7QcnzxxRe6/fbb9fXXX+uSSy7RDTfcoF27dumSSy6RJD377LMKDQ1VVlaW3G63MjIy9OKLLwZ51AAABA5BGwAgSUrs9O31sc9NGqrTXQcZ6bNzxwhdGt+0m2Oh5Vq9evV5t0dGRqqgoEAFBQUBGlHbkKJjijy2Vwpp/hkqkqToLlI8N5QDgGAgaAMAvPS+pJOU0vwj5AAuXoeqL7XF9pCi17ov3PiiO42Wct4lbANAEBC0AQAAgizs9HFFh7h1aORzSr1iaPM7PPYPac006eTXBG0ACAKCNgAAQAvhju8tpQwN9jAAAM0UGuwBAAAAAADQlhC0AQAAAAAwiKANAAAAAIBBBG0AAAAAAAziZmgA0IoZfe7usX80vw8AAAAQtAGgtfLbc3eju5jrDwAAoB0iaANAK2X8ubvStyGbZ+4CF8/QmSC2ilIj/QAAWgaCNgC0cjx3FwiC6C7fngGyZpqR7lIlnbRsqo9MMNIfACC4CNoAAAC+ik+Vct6VTn5tpLvSr6p01x8P6PedLjXSHwAguAjaAAAATRGfauxSi9NWpQ6r0khfAIDg4/FeAAAAAAAYRNAGAAAAAMAggjYAAAAAAAYRtAEAAAAAMIigDQAAAACAQcaDdn5+vq655hrFxMQoMTFRt956q0pKSrzajBgxQiEhIV7LfffdZ3ooAAAAAAAEnPHHe23fvl05OTm65pprVFdXp1//+tcaM2aM9u/fr44dO3raTZs2TfPmzfO8jo6ONj0UAACAdq30qyqdtsw8NqxzxwhdGh9lpC8AaOuMB+2NGzd6vV6+fLkSExNVVFSkm266ybM+Ojpadrv9ovp0u91yu92e1y6Xy8xgAQAA2qCjVW4lSvr56mLtMxS0ozqEacsvbiZsA8BFMB60v6+y8ttf7gkJCV7rV65cqVdeeUV2u10TJkzQI488cs6j2vn5+Zo7d66/hwoAANAmuE7VKlHSL8f00SVXXtvs/kqPVunBV4v1TXUNQRsALoJfg3ZDQ4MefPBBXX/99Ro4cKBn/R133KHu3bsrJSVFH374oWbNmqWSkhKtWbPmrP3k5eUpNzfX89rlcik1NdWfQwcAAAi40qNVRvr56vgp9ZaUmhCl3pfGGekTAHDx/Bq0c3Jy9NFHH+ntt9/2Wj99+nTPvwcNGqTk5GSNHj1aBw4cUK9evc7ox2azyWaz+XOoAAAAQdO5Y4SiOoTpwVeLjfQ3IOSgRtqk2KgORvoDAPjGb0F7xowZWr9+vXbs2KHLLrvsvG3T0tIkSaWlpWcN2gDQFnxZcUrfVNcY66/xiBUQaPn5+VqzZo3+/ve/KyoqSj/4wQ/029/+Vn369PG0GTFihLZv3+71vnvvvVdLliwJ9HBbhUvjo7TlFzcb+x0ReSxOWislduJABQAEg/GgbVmWHnjgAa1du1bbtm1Tjx49Lvie4uJiSVJycrLp4QBAi/BlxSlN+d3/p6i6CmN99g75UiMjOGKFwOMJI/5xaXyUueufQzqZ6QcA0CTGg3ZOTo5WrVqlP//5z4qJiZHT6ZQkxcXFKSoqSgcOHNCqVas0btw4denSRR9++KFmzpypm266SYMHDzY9HABoEarKD2p9aK6ibe4LN/ZBQ3iUEhNTjPYJXIg/njACPzn2DyPdRB6rUoqOGekLANoD40F78eLFkr49Zey7li1bprvvvlsRERHasmWLFi1apOrqaqWmpiorK0uzZ882PRQAaDHCTh9XdIhbh0Y+p9QrhhrrNzS6ixTPzSERXCaeMMKjPA2L7iJ1iJbWTDPSXW9JW2w2Haq6RhI3VwOAC/HLqePnk5qaesY1WwDQXrjje0spQ4M9DMAYU08Y4VGehsWnSjnvSie/NtLdoU+KlfrWzxV2+riR/gCgrfP7c7QBAEDbZeoJIzzK0w/iU42d8eL+ysxjxwCgvSBoA7h4FYeMHR2R9O2pjS34tGeTdwnnDuFoi0w+YYRHeQIA2hKCNoCLU3FIKrhWqj1prs8O0d+e2tgCw/aXFaeU/rvtOlVbb6Q/nmmLtoQnjAAAcH4EbQAX5+TX34bsiX+Qul7Z/P6O/ePbm/Sc/LpFBu1vqmt0qrZei24bqt6JzX9MDs+0RVvCE0YAADg/gjYA33S9sl3dzKt3YicNvNTAHXZ5pi3aEJ4w0n7ZKkqlwwZ/n7XwS4gAoKkI2pL5604lCgcAoM3iCSPtT31kgk5aNqW+9XPpLXP9NoRHKXTGHv5mAtDmELT9cd2p1KKvPYUh7ezGYH5z7B/m+mqvcwgAftYpqYf+T8MziqqrMNZn75Av9Zxe1NGjh5XI724AbQxB2/R1p1KLv/YUBrSzG4P5RXSXb/d5zTRzfba3OQSAALk0Pkr/9xdZxp7EIElf/eNdafuLcp2qVaKxXtFsHEgAjCBoN2pn152imfx1Y7CyQnPFraUXtvjUb0Oxqf01PIeRx6qUomMGBgYAbcOl8VG6ND7KWH+lx8z15TftLXRyIAEwhqANNIepL2ja69Hd+FRz4zM8h70lbbHZdKjqGkkGboYGAGhd2mPobGdPGAH8iaANtAT+OrrbngpbfKqcU3ao6ptyI91VfP6Rhr8/S2GnjxvpDwDQyrTn0MmZnkCzEbSBlsLk0d126MuKU0r/wwGdqq030t+AkDq9bpNiozoY6Q8A0EoROoHAawOXbRC0AbQJ31TX6FRtvRbdNlS9E5v/jNfIY3HSWimxk83A6AAAAWHyj3OTT8UAcPHayGUbBG2gLTP5R0Ir+YOjd2InDbzUwDXVIf8vrJva71YyfwDQavnrj/PoLub6A3BhbeSyDYI20Bb54+ZqUvv6g8NfN6hrL/MHABfp0PFTOv1lZbP7iTz2uXqbfmRrS79LONCWtfLLNgja7VUbuO4B52H65mqN2tP/z/6Yw/Y0fwBwAY33wHj6jRLt29T853MPCDmo121S0clLZLN6NLs/SeqsCF1qpCcA7Q1Buz1qI9c94AK4uVrzMYcA4DeN98B4btJQne46qNn9uQ9FSxulOX/ep32Wmb9xojqEacsvbjb6/HDjlxLxJS7QIhG026M2ct0DAABo/Xpf0klKMXFvjc6SzAX30qNVevDVYn1TXWMmaPvzsi4OdgAtDkG7PWvl1z0AAAB8n7Hgbpo/LklqPNhRVmimX27c2TJxyWerRNAGAABA8LSnpzuYviSJG3e2fVzy2WoRtP3J5C98vnkCAABtCSGx+bhxZ9vHJZ+tFkHbH/xVOPjmCQAAtBWERDNaw407OfjUfFzy2eoQtP3BdOHg+hsAANAWtYKQWHq0ylhfnTtGmL2DeUvHwaf2wfQ15G0kqxC0/cVk4WiPp1aZ/g/WtDbyCyDYvqw4pW+qm//sVMnsH0IAAHTuGKGoDmF68NViY3365XFhLZm/Dj5x2nPL4Y9ryKWWn1UuAkG7NWhvp1b56z9Y09rAL4Bg+rLilNJ/t12nauuN9RnVIUydO0YY6w8A0H5dGh+lLb+42egXwkYfF9ZatIKzFtAMpq8hb9SSs8pFImi3Fq3hl5TJu4b64z9Y09rAL4Bg+qa6Rqdq67XotqHqndjJSJ/t7pQ8AIBfXRofRV2B70yememvsyhN3+2fa8jPQNBG8/nr1PZuDoJsO9A7sZMGXtoCn3cKAADaB5PB+OQx6dUp5h/HZeosyvZ4SWqQELTRfO3t1PZWwuT1z/7ANdUAACDo/PWc6jv/Pym6q5n+TP5dzN/tAUPQhhmt4dT2Fs5kMP66ukb3/d8io9c/+wPXVAMA0Dymv1hvFZdhmTyd2h+XLLb04Mnf7QFB0AZaAH/dGGzFPdeqSwsOsq2imAMA0EL56++HJVOGGfv7wWit98dpzxKXLMIvghq0CwoK9NRTT8npdGrIkCF64YUXdO211wZzSMBFM/1oKm4MBqAtotYD/mP6xqKNZ8Rlv/SugdF9y+gjzfxx2rPU8o9Ao1UKWtB+9dVXlZubqyVLligtLU2LFi1SRkaGSkpKlJiYGKxhARfFX98gX9MjgXAMoM2g1gOBYfLGoi3+kWac9oxWImhB+5lnntG0adP0k5/8RJK0ZMkSvf7663rppZf0q1/9KljDAi4Kj6YCgAuj1gNnMnkzUH/cWNQfjzQzOU7+XkJrEZSgXVNTo6KiIuXl5XnWhYaGKj09XYWFhWe0d7vdcrvdnteVlZWSJJfL1fzBnKiS3JYOHCyX23Wo+f2hXfj0q2o1uE/KHtWgbjEhhnqtlctVa6gvABfjRFW1XG7r2/9tZk1pfL9lWSaG1uq1qFovs/9fA00RXn9aEQ2n9bOXdxrtN7JDqMLrT8vlMvX3iDn+2OfIDqFaNOkqJUR3MNYn2jbb1+Xq5ba+zX0BrPVBCdrHjh1TfX29kpKSvNYnJSXp73//+xnt8/PzNXfu3DPWp6YaPG1kwVhzfaHdcCwK9ggAGGGwBpw4cUJxcTwbvkXWeol6jzap31PBHkFgjW9n+wtDFtxorKuLqfWt4q7jeXl5ys3N9bxuaGjQ8ePH1aVLF4WENO/bO5fLpdTUVB06dEixsbHNHWq7wJz5jjnzHXPmO+bMdybnzLIsnThxQikpKYZG1774s9ZL/PfhK+bLd8yZ75gz3zFnvgtWrQ9K0O7atavCwsJUXl7utb68vFx2u/2M9jabTTabzWtdfHy80THFxsbyw+oj5sx3zJnvmDPfMWe+MzVnHMn+p5ZY6yX++/AV8+U75sx3zJnvmDPfBbrWhzb7k5ogIiJCw4YN09atWz3rGhoatHXrVjkcjmAMCQAAGEStBwC0Z0E7dTw3N1fZ2dkaPny4rr32Wi1atEjV1dWeO5MCAIDWjVoPAGivgha0b7vtNn311VeaM2eOnE6nhg4dqo0bN55x0xR/s9lsevTRR884XQ3nxpz5jjnzHXPmO+bMd8yZf7WUWi/x/7WvmC/fMWe+Y858x5z5LlhzFmLxHBIAAAAAAIwJyjXaAAAAAAC0VQRtAAAAAAAMImgDAAAAAGAQQRsAAAAAAIMI2gAAAAAAGNQugnZBQYEuv/xyRUZGKi0tTe++++5527/22mvq27evIiMjNWjQIG3YsCFAI205fJmzP/zhD7rxxhvVuXNnde7cWenp6Rec47bI15+zRqtXr1ZISIhuvfVW/w6wBfJ1zioqKpSTk6Pk5GTZbDZdeeWV7e6/T1/nbNGiRerTp4+ioqKUmpqqmTNn6vTp0wEabXDt2LFDEyZMUEpKikJCQrRu3boLvmfbtm26+uqrZbPZ1Lt3by1fvtzv44Q51HvfUOt9R633HbXed9T6i9eia73Vxq1evdqKiIiwXnrpJWvfvn3WtGnTrPj4eKu8vPys7d955x0rLCzMWrhwobV//35r9uzZVocOHay9e/cGeOTB4+uc3XHHHVZBQYH1wQcfWB9//LF19913W3FxcdYXX3wR4JEHj69z1ujgwYPWpZdeat14443WLbfcEpjBthC+zpnb7baGDx9ujRs3znr77betgwcPWtu2bbOKi4sDPPLg8XXOVq5cadlsNmvlypXWwYMHrU2bNlnJycnWzJkzAzzy4NiwYYP1m9/8xlqzZo0lyVq7du1523/66adWdHS0lZuba+3fv9964YUXrLCwMGvjxo2BGTCahXrvG2q976j1vqPW+45a75uWXOvbfNC+9tprrZycHM/r+vp6KyUlxcrPzz9r+x//+MfW+PHjvdalpaVZ9957r1/H2ZL4OmffV1dXZ8XExFgrVqzw1xBbnKbMWV1dnfWDH/zA+q//+i8rOzu73RVfX+ds8eLFVs+ePa2amppADbHF8XXOcnJyrFGjRnmty83Nta6//nq/jrMlupji+/DDD1sDBgzwWnfbbbdZGRkZfhwZTKHe+4Za7ztqve+o9b6j1jddS6v1bfrU8ZqaGhUVFSk9Pd2zLjQ0VOnp6SosLDzrewoLC73aS1JGRsY527c1TZmz7zt58qRqa2uVkJDgr2G2KE2ds3nz5ikxMVFTp04NxDBblKbM2V/+8hc5HA7l5OQoKSlJAwcO1JNPPqn6+vpADTuomjJnP/jBD1RUVOQ55ezTTz/Vhg0bNG7cuICMubVp77//WzPqvW+o9b6j1vuOWu87ar3/BfJ3f7jxHluQY8eOqb6+XklJSV7rk5KS9Pe///2s73E6nWdt73Q6/TbOlqQpc/Z9s2bNUkpKyhk/xG1VU+bs7bff1tKlS1VcXByAEbY8TZmzTz/9VG+++aYmT56sDRs2qLS0VD/96U9VW1urRx99NBDDDqqmzNkdd9yhY8eO6YYbbpBlWaqrq9N9992nX//614EYcqtzrt//LpdLp06dUlRUVJBGhguh3vuGWu87ar3vqPW+o9b7XyBrfZs+oo3AW7BggVavXq21a9cqMjIy2MNpkU6cOKEpU6boD3/4g7p27Rrs4bQaDQ0NSkxM1O9//3sNGzZMt912m37zm99oyZIlwR5ai7Vt2zY9+eSTevHFF/X+++9rzZo1ev311zV//vxgDw1AK0atvzBqfdNQ631HrW+52vQR7a5duyosLEzl5eVe68vLy2W328/6Hrvd7lP7tqYpc9bo6aef1oIFC7RlyxYNHjzYn8NsUXydswMHDuizzz7ThAkTPOsaGhokSeHh4SopKVGvXr38O+gga8rPWXJysjp06KCwsDDPun79+snpdKqmpkYRERF+HXOwNWXOHnnkEU2ZMkX//u//LkkaNGiQqqurNX36dP3mN79RaCjftX7XuX7/x8bGcjS7haPe+4Za7ztqve+o9b6j1vtfIGt9m575iIgIDRs2TFu3bvWsa2ho0NatW+VwOM76HofD4dVekjZv3nzO9m1NU+ZMkhYuXKj58+dr48aNGj58eCCG2mL4Omd9+/bV3r17VVxc7Fl+9KMfaeTIkSouLlZqamoghx8UTfk5u/7661VaWur5Q0WS/vGPfyg5ObnNF16paXN28uTJMwps4x8v394zBN/V3n//t2bUe99Q631Hrfcdtd531Hr/C+jvfuO3V2thVq9ebdlsNmv58uXW/v37renTp1vx8fGW0+m0LMuypkyZYv3qV7/ytH/nnXes8PBw6+mnn7Y+/vhj69FHH21Xj/uwLN/nbMGCBVZERIT13//939aRI0c8y4kTJ4K1CwHn65x9X3u8E6mvc1ZWVmbFxMRYM2bMsEpKSqz169dbiYmJ1uOPPx6sXQg4X+fs0UcftWJiYqw//vGP1qeffmq98cYbVq9evawf//jHwdqFgDpx4oT1wQcfWB988IElyXrmmWesDz74wPr8888ty7KsX/3qV9aUKVM87Rsf+fHQQw9ZH3/8sVVQUMDjvVoR6r1vqPW+o9b7jlrvO2q9b1pyrW/zQduyLOuFF16wunXrZkVERFjXXnuttWvXLs+2m2++2crOzvZq/6c//cm68sorrYiICGvAgAHW66+/HuARB58vc9a9e3dL0hnLo48+GviBB5GvP2ff1R6Lr2X5Pmc7d+600tLSLJvNZvXs2dN64oknrLq6ugCPOrh8mbPa2lrrscces3r16mVFRkZaqamp1k9/+lPrm2++CfzAg+Ctt9466++mxjnKzs62br755jPeM3ToUCsiIsLq2bOntWzZsoCPG01HvfcNtd531HrfUet9R62/eC251odYFucUAAAAAABgSpu+RhsAAAAAgEAjaAMAAAAAYBBBGwAAAAAAgwjaAAAAAAAYRNAGAAAAAMAggjYAAAAAAAYRtAEAAAAAMIigDQAAAACAQQRtAAAAAAAMImgDAAAAAGAQQRsAAAAAAIMI2gAAAAAAGETQBgAAAADAIII2AAAAAAAGEbQBAAAAADCIoA0AAAAAgEEEbQAAAAAADCJoAwAAAABgEEEbAAAAAACDCNoAAAAAABhE0AYAAAAAwCCCNgAAAAAABhG0AQAAAAAwiKANAAAAAIBBBG0AAAAAAAwiaAMAAAAAYBBBGwAAAAAAgwjaAAAAAAAYRNAGAAAAAMAggjYAAAAAAAYRtAEAAAAAMIigDQAAAACAQQRtAAAAAAAMImgDAAAAAGAQQRsAAAAAAIMI2gAAAAAAGETQBgAAAADAIII2AAAAAAAGEbQBAAAAADAoPNgDaIqGhgYdPnxYMTExCgkJCfZwAACQZVk6ceKEUlJSFBrK99jNRa0HALQ0vtT6Vhm0Dx8+rNTU1GAPAwCAMxw6dEiXXXZZsIfR6lHrAQAt1cXU+lYZtGNiYiR9u4OxsbFBHg0AAJLL5VJqaqqnRrUVCxYsUF5enn7+859r0aJFkqTTp0/rF7/4hVavXi23262MjAy9+OKLSkpK8ryvrKxM999/v9566y116tRJ2dnZys/PV3j4xf3pQa0HALQ0vtT6ZgXtYBXfxlPIYmNjKb4AgBalLZ3mvGfPHv3nf/6nBg8e7LV+5syZev311/Xaa68pLi5OM2bM0MSJE/XOO+9Ikurr6zV+/HjZ7Xbt3LlTR44c0V133aUOHTroySefvKjPptYDAFqqi6n1Tb6I7HzF969//atee+01bd++XYcPH9bEiRM92xuLb01NjXbu3KkVK1Zo+fLlmjNnTlOHAgAADKuqqtLkyZP1hz/8QZ07d/asr6ys1NKlS/XMM89o1KhRGjZsmJYtW6adO3dq165dkqQ33nhD+/fv1yuvvKKhQ4cqMzNT8+fPV0FBgWpqas76eW63Wy6Xy2sBAKC1alLQDnTxBQAAgZWTk6Px48crPT3da31RUZFqa2u91vft21fdunVTYWGhJKmwsFCDBg3yOpstIyNDLpdL+/btO+vn5efnKy4uzrNwfTYAoDVrUtAOdPHlW24AAAJn9erVev/995Wfn3/GNqfTqYiICMXHx3utT0pKktPp9LT5bp1v3N647Wzy8vJUWVnpWQ4dOmRgTwAACA6fr9FuLL579uw5Y5u/im9+fr7mzp3r61ABAICPDh06pJ///OfavHmzIiMjA/a5NptNNpstYJ8HAIA/+XREu7H4rly5MqDFl2+5AQAIjKKiIh09elRXX321wsPDFR4eru3bt+v5559XeHi4kpKSVFNTo4qKCq/3lZeXy263S5LsdrvKy8vP2N64DQCAts6noB2s4muz2Tx3HeXuowAA+M/o0aO1d+9eFRcXe5bhw4dr8uTJnn936NBBW7du9bynpKREZWVlcjgckiSHw6G9e/fq6NGjnjabN29WbGys+vfvH/B9AgAg0Hw6dbyx+H7XT37yE/Xt21ezZs1Samqqp/hmZWVJOnvxfeKJJ3T06FElJiZKovgCANBSxMTEaODAgV7rOnbsqC5dunjWT506Vbm5uUpISFBsbKweeOABORwOXXfddZKkMWPGqH///poyZYoWLlwop9Op2bNnKycnh9PDAQDtgk9Bm+ILAACeffZZhYaGKisrS263WxkZGXrxxRc928PCwrR+/Xrdf//9cjgc6tixo7KzszVv3rwgjhoAgMDx+WZoF0LxBQCgbdm2bZvX68jISBUUFKigoOCc7+nevbs2bNjg55EBANAyhViWZQV7EL5yuVyKi4tTZWUl12sDAFoEapNZzCcAoKXxpTY16TnaAAAAAADg7IyfOt4afVlxSt9U1xjts3PHCF0aH2W0TwAA0AwVh6STX5vrL7qLFJ9qrj8AQJvR7oP2lxWnlP677TpVW2+036gOYdryi5sJ2wAAtAQVh6SCa6Xak+b67BAt5bxL2AYAnKHdB+1vqmt0qrZei24bqt6JnYz0WXq0Sg++WqxvqmsI2gAAtAQnv/42ZE/8g9T1yub3d+wf0ppp3/ZL0AYAfE+7D9qNeid20sBL44I9DAAA4EelVopOWz2a3U+kVaXeBsYDAGibCNoAAKDNO1rlVqKkn68u1j6rstn9DQg5qNdt/+wXAIDvImgDAIA2z3WqVomSfjmmjy658tpm9/fVPyKk7f/sFwCA7yJoS0rRMUUe2yuFmLlGO/JYlVJ0zEhfAADAnNSEKPU2cKlY6THuwQIAOLd2H7Q7VH2pLbaHFL3WbazP3pK22Gw6VHWNJK77BgAAAID2pN0H7bDTxxUd4tahkc8p9YqhRvo89EmxUt/6ucJOHzfSHwAAAACg9Wj3QbuRO763lDLUTF9fVRnpBwAAAADQ+oQGewAAAAAAALQlBG0AAAAAAAwiaAMAAAAAYBBBGwAAAAAAgwjaAAAAAAAYRNAGAAAAAMAggjYAAAAAAAYRtAEAAAAAMIigDQAAAACAQQRtAAAAAAAMImgDAAAAAGAQQRsAAAAAAIMI2gAAAAAAGETQBgAAAADAIII2AAAAAAAGEbQBAIDH4sWLNXjwYMXGxio2NlYOh0P/8z//49k+YsQIhYSEeC333XefVx9lZWUaP368oqOjlZiYqIceekh1dXWB3hUAAIImPNgDAAAALcdll12mBQsW6IorrpBlWVqxYoVuueUWffDBBxowYIAkadq0aZo3b57nPdHR0Z5/19fXa/z48bLb7dq5c6eOHDmiu+66Sx06dNCTTz4Z8P0BACAYfDqizbfcAAC0bRMmTNC4ceN0xRVX6Morr9QTTzyhTp06adeuXZ420dHRstvtniU2Ntaz7Y033tD+/fv1yiuvaOjQocrMzNT8+fNVUFCgmpqaYOwSAAAB51PQbvyWu6ioSO+9955GjRqlW265Rfv27fO0mTZtmo4cOeJZFi5c6NnW+C13TU2Ndu7cqRUrVmj58uWaM2eOuT0CAABG1NfXa/Xq1aqurpbD4fCsX7lypbp27aqBAwcqLy9PJ0+e9GwrLCzUoEGDlJSU5FmXkZEhl8vl9ffC97ndbrlcLq8FAIDWyqdTxydMmOD1+oknntDixYu1a9cuz+lkjd9yn03jt9xbtmxRUlKShg4dqvnz52vWrFl67LHHFBER0cTdAAAApuzdu1cOh0OnT59Wp06dtHbtWvXv31+SdMcdd6h79+5KSUnRhx9+qFmzZqmkpERr1qyRJDmdTq+QLcnz2ul0nvMz8/PzNXfuXD/tEQAAgdXkm6HxLTcAAG1Tnz59VFxcrN27d+v+++9Xdna29u/fL0maPn26MjIyNGjQIE2ePFkvv/yy1q5dqwMHDjTrM/Py8lRZWelZDh06ZGJXAAAICp9vhsa33AAAtG0RERHq3bu3JGnYsGHas2ePnnvuOf3nf/7nGW3T0tIkSaWlperVq5fsdrveffddrzbl5eWSdM4z3iTJZrPJZrOZ2gUAAILK56Dd+C13ZWWl/vu//1vZ2dnavn27+vfvr+nTp3vaDRo0SMnJyRo9erQOHDigXr16NXmQeXl5ys3N9bx2uVxKTU1tcn8AAODiNTQ0yO12n3VbcXGxJCk5OVmS5HA49MQTT+jo0aNKTEyUJG3evFmxsbGeL+YBAGjrfA7afMsNAEDblZeXp8zMTHXr1k0nTpzQqlWrtG3bNm3atEkHDhzQqlWrNG7cOHXp0kUffvihZs6cqZtuukmDBw+WJI0ZM0b9+/fXlClTtHDhQjmdTs2ePVs5OTnUcgBAu9Hka7Qb+fot9969e3X06FFPG77lBgCg5Th69Kjuuusu9enTR6NHj9aePXu0adMm/fCHP1RERIS2bNmiMWPGqG/fvvrFL36hrKws/fWvf/W8PywsTOvXr1dYWJgcDofuvPNO3XXXXV7P3QYAoK3z6Yg233IDANC2LV269JzbUlNTtX379gv20b17d23YsMHksAAAaFV8CtqN33IfOXJEcXFxGjx4sOdb7kOHDmnLli1atGiRqqurlZqaqqysLM2ePdvz/sZvue+//345HA517NhR2dnZfMsNAAAAAGgzfArafMsNAAAAAMD5NfsabQAAAAAA8E8EbQAAAAAADCJoAwAAAABgEEEbAAAAAACDCNoAAAAAABhE0AYAAAAAwCCCNgAAAAAABhG0AQAAAAAwiKANAAAAAIBBBG0AAAAAAAwiaAMAAAAAYBBBGwAAAAAAgwjaAAAAAAAYRNAGAAAAAMAggjYAAAAAAAYRtAEAAAAAMIigDQAAAACAQQRtAAAAAAAMImgDAAAAAGAQQRsAAAAAAIMI2gAAAAAAGETQBgAAAADAIII2AAAAAAAGEbQBAAAAADCIoA0AAAAAgEEEbQAAAAAADCJoAwAAj8WLF2vw4MGKjY1VbGysHA6H/ud//sez/fTp08rJyVGXLl3UqVMnZWVlqby83KuPsrIyjR8/XtHR0UpMTNRDDz2kurq6QO8KAABBQ9AGAAAel112mRYsWKCioiK99957GjVqlG655Rbt27dPkjRz5kz99a9/1Wuvvabt27fr8OHDmjhxouf99fX1Gj9+vGpqarRz506tWLFCy5cv15w5c4K1SwAABJxPQZtvuQEAaNsmTJigcePG6YorrtCVV16pJ554Qp06ddKuXbtUWVmppUuX6plnntGoUaM0bNgwLVu2TDt37tSuXbskSW+88Yb279+vV155RUOHDlVmZqbmz5+vgoIC1dTUBHnvAAAIDJ+CNt9yAwDQftTX12v16tWqrq6Ww+FQUVGRamtrlZ6e7mnTt29fdevWTYWFhZKkwsJCDRo0SElJSZ42GRkZcrlcnr8XzsbtdsvlcnktAAC0VuG+NJ4wYYLX6yeeeEKLFy/Wrl27dNlll2np0qVatWqVRo0aJUlatmyZ+vXrp127dum6667zfMu9ZcsWJSUlaejQoZo/f75mzZqlxx57TBEREWf9XLfbLbfb7XlN8QUAwH/27t0rh8Oh06dPq1OnTlq7dq369++v4uJiRUREKD4+3qt9UlKSnE6nJMnpdHqF7MbtjdvOJT8/X3PnzjW7IwAABEmTr9EO5Lfc+fn5iouL8yypqalNHTYAALiAPn36qLi4WLt379b999+v7Oxs7d+/36+fmZeXp8rKSs9y6NAhv34eAAD+5HPQ3rt3rzp16iSbzab77rvP8y230+n027fcFF8AAAInIiJCvXv31rBhw5Sfn68hQ4boueeek91uV01NjSoqKrzal5eXy263S5LsdvsZ92dpfN3Y5mxsNpvnHjCNCwAArZXPQTsY33JTfAEACJ6Ghga53W4NGzZMHTp00NatWz3bSkpKVFZWJofDIUlyOBzau3evjh496mmzefNmxcbGqn///gEfOwAAweDTNdrSP7/llqRhw4Zpz549eu6553Tbbbd5vuX+7lHt73/L/e6773r1dzHfcgMAgMDIy8tTZmamunXrphMnTmjVqlXatm2bNm3apLi4OE2dOlW5ublKSEhQbGysHnjgATkcDl133XWSpDFjxqh///6aMmWKFi5cKKfTqdmzZysnJ0c2my3IewcAQGA0+znafMsNAEDbcfToUd11113q06ePRo8erT179mjTpk364Q9/KEl69tln9X/+z/9RVlaWbrrpJtntdq1Zs8bz/rCwMK1fv15hYWFyOBy68847ddddd2nevHnB2iUAAALOpyPafMsNAEDbtnTp0vNuj4yMVEFBgQoKCs7Zpnv37tqwYYPpoQEA0Gr4FLQbv+U+cuSI4uLiNHjw4DO+5Q4NDVVWVpbcbrcyMjL04osvet7f+C33/fffL4fDoY4dOyo7O5tvuQEAAAAAbYZPQZtvuQEAAAAAOL9mX6MNAAAAAAD+iaANAAAAAIBBBG0AAAAAAAwiaAMAAAAAYBBBGwAAAAAAgwjaAAAAAAAYRNAGAAAAAMAggjYAAAAAAAYRtAEAAAAAMIigDQAAAACAQQRtAAAAAAAMImgDAAAAAGAQQRsAAAAAAIMI2gAAAAAAGETQBgAAAADAIII2AAAAAAAGEbQBAAAAADCIoA0AAAAAgEEEbQAAAAAADCJoAwAAAABgEEEbAAAAAACDCNoAAAAAABhE0AYAAAAAwCCCNgAAAAAABhG0AQAAAAAwiKANAAAAAIBBBG0AAOCRn5+va665RjExMUpMTNStt96qkpISrzYjRoxQSEiI13Lfffd5tSkrK9P48eMVHR2txMREPfTQQ6qrqwvkrgAAEDQ+BW2KLwAAbdv27duVk5OjXbt2afPmzaqtrdWYMWNUXV3t1W7atGk6cuSIZ1m4cKFnW319vcaPH6+amhrt3LlTK1as0PLlyzVnzpxA7w4AAEER7kvjxuJ7zTXXqK6uTr/+9a81ZswY7d+/Xx07dvS0mzZtmubNm+d5HR0d7fl3Y/G12+3auXOnjhw5orvuuksdOnTQk08+aWCXAABAU23cuNHr9fLly5WYmKiioiLddNNNnvXR0dGy2+1n7eONN97Q/v37tWXLFiUlJWno0KGaP3++Zs2apccee0wRERF+3QcAAILNpyPaGzdu1N13360BAwZoyJAhWr58ucrKylRUVOTVrrH4Ni6xsbGebY3F95VXXtHQoUOVmZmp+fPnq6CgQDU1NWb2CgAAGFFZWSlJSkhI8Fq/cuVKde3aVQMHDlReXp5Onjzp2VZYWKhBgwYpKSnJsy4jI0Mul0v79u076+e43W65XC6vBQCA1qpZ12hTfAEAaLsaGhr04IMP6vrrr9fAgQM96++44w698soreuutt5SXl6f/+3//r+68807PdqfT6VXnJXleO53Os35Wfn6+4uLiPEtqaqof9ggAgMDw6dTx7zpf8e3evbtSUlL04YcfatasWSopKdGaNWskNb34zp07t6lDBQAATZCTk6OPPvpIb7/9ttf66dOne/49aNAgJScna/To0Tpw4IB69erVpM/Ky8tTbm6u57XL5SJsAwBarSYHbYovAABt14wZM7R+/Xrt2LFDl1122XnbpqWlSZJKS0vVq1cv2e12vfvuu15tysvLJemc13XbbDbZbDYDIwcAIPiadOp4Y/F96623fCq+0rcFtrHYNrqY4hsbG+u1AAAA8yzL0owZM7R27Vq9+eab6tGjxwXfU1xcLElKTk6WJDkcDu3du1dHjx71tNm8ebNiY2PVv39/v4wbAICWxKegTfEFAKBty8nJ0SuvvKJVq1YpJiZGTqdTTqdTp06dkiQdOHBA8+fPV1FRkT777DP95S9/0V133aWbbrpJgwcPliSNGTNG/fv315QpU/S///u/2rRpk2bPnq2cnByOWgMA2gWfTh3PycnRqlWr9Oc//9lTfCUpLi5OUVFROnDggFatWqVx48apS5cu+vDDDzVz5sxzFt+FCxfK6XRSfAEAaCEWL14sSRoxYoTX+mXLlunuu+9WRESEtmzZokWLFqm6ulqpqanKysrS7NmzPW3DwsK0fv163X///XI4HOrYsaOys7O9Hv0JAEBb5lPQpvgCANC2WZZ13u2pqanavn37Bfvp3r27NmzYYGpYAAC0Kj4FbYovAAAAAADn16znaAMAAAAAAG8EbQAAAAAADCJoAwAAAABgEEEbAAAAAACDCNoAAAAAABhE0AYAAAAAwCCCNgAAAAAABhG0AQAAAAAwiKANAAAAAIBBBG0AAAAAAAwiaAMAAAAAYBBBGwAAAAAAgwjaAAAAAAAYRNAGAAAAAMAggjYAAAAAAAYRtAEAAAAAMIigDQAAAACAQQRtAAAAAAAMImgDAAAAAGAQQRsAAAAAAIMI2gAAAAAAGETQBgAAAADAIII2AAAAAAAGEbQBAAAAADCIoA0AAAAAgEEEbQAAAAAADCJoAwAAAABgEEEbAAB45Ofn65prrlFMTIwSExN16623qqSkxKvN6dOnlZOToy5duqhTp07KyspSeXm5V5uysjKNHz9e0dHRSkxM1EMPPaS6urpA7goAAEHjU9Cm+AIA0LZt375dOTk52rVrlzZv3qza2lqNGTNG1dXVnjYzZ87UX//6V7322mvavn27Dh8+rIkTJ3q219fXa/z48aqpqdHOnTu1YsUKLV++XHPmzAnGLgEAEHA+BW2KLwAAbdvGjRt19913a8CAARoyZIiWL1+usrIyFRUVSZIqKyu1dOlSPfPMMxo1apSGDRumZcuWaefOndq1a5ck6Y033tD+/fv1yiuvaOjQocrMzNT8+fNVUFCgmpqaYO4eAAAB4VPQpvgCANC+VFZWSpISEhIkSUVFRaqtrVV6erqnTd++fdWtWzcVFhZKkgoLCzVo0CAlJSV52mRkZMjlcmnfvn1n/Ry32y2Xy+W1AADQWjXrGm2KLwAAbVdDQ4MefPBBXX/99Ro4cKAkyel0KiIiQvHx8V5tk5KS5HQ6PW2+W+cbtzduO5v8/HzFxcV5ltTUVMN7AwBA4DQ5aFN8AQBo23JycvTRRx9p9erVfv+svLw8VVZWepZDhw75/TMBAPCXJgdtii8AAG3XjBkztH79er311lu67LLLPOvtdrtqampUUVHh1b68vFx2u93T5vs3Qm183djm+2w2m2JjY70WAABaqyYFbYovAABtk2VZmjFjhtauXas333xTPXr08No+bNgwdejQQVu3bvWsKykpUVlZmRwOhyTJ4XBo7969Onr0qKfN5s2bFRsbq/79+wdmRwAACCKfgjbFFwCAti0nJ0evvPKKVq1apZiYGDmdTjmdTp06dUqSFBcXp6lTpyo3N1dvvfWWioqK9JOf/EQOh0PXXXedJGnMmDHq37+/pkyZov/93//Vpk2bNHv2bOXk5MhmswVz9wAACIhwXxrn5ORo1apV+vOf/+wpvtK3RTcqKsqr+CYkJCg2NlYPPPDAOYvvwoUL5XQ6Kb4AALQQixcvliSNGDHCa/2yZct09913S5KeffZZhYaGKisrS263WxkZGXrxxRc9bcPCwrR+/Xrdf//9cjgc6tixo7KzszVv3rxA7QYAAEHlU9Cm+AIA0LZZlnXBNpGRkSooKFBBQcE523Tv3l0bNmwwOTQAAFoNn4I2xRcAAAAAgPNr1nO0AQAAAACAN4I2AAAAAAAGEbQBAAAAADCIoA0AAAAAgEEEbQAAAAAADCJoAwAAAABgEEEbAAAAAACDCNoAAAAAABhE0AYAAAAAwCCCNgAAAAAABhG0AQAAAAAwiKANAAAAAIBBBG0AAAAAAAwiaAMAAAAAYBBBGwAAAAAAgwjaAAAAAAAYRNAGAAAAAMAggjYAAAAAAAYRtAEAAAAAMIigDQAAAACAQQRtAAAAAAAMImgDAAAAAGAQQRsAAAAAAIMI2gAAAAAAGETQBgAAAADAIII2AAAAAAAGEbQBAAAAADCIoA0AALzs2LFDEyZMUEpKikJCQrRu3Tqv7XfffbdCQkK8lrFjx3q1OX78uCZPnqzY2FjFx8dr6tSpqqqqCuBeAAAQPD4HbYovAABtW3V1tYYMGaKCgoJzthk7dqyOHDniWf74xz96bZ88ebL27dunzZs3a/369dqxY4emT5/u76EDANAihPv6hsbie88992jixIlnbTN27FgtW7bM89pms3ltnzx5so4cOaLNmzertrZWP/nJTzR9+nStWrXK1+EAAADDMjMzlZmZed42NptNdrv9rNs+/vhjbdy4UXv27NHw4cMlSS+88ILGjRunp59+WikpKcbHDABAS+Jz0A5G8XW73XK73Z7XLpfL12EDAACDtm3bpsTERHXu3FmjRo3S448/ri5dukiSCgsLFR8f76nzkpSenq7Q0FDt3r1b//Iv/3JGf9R6AEBb4pdrtBuLb58+fXT//ffr66+/9my7UPE9m/z8fMXFxXmW1NRUfwwbAABchLFjx+rll1/W1q1b9dvf/lbbt29XZmam6uvrJUlOp1OJiYle7wkPD1dCQoKcTudZ+6TWAwDaEuNB2x/FNy8vT5WVlZ7l0KFDpocNAAAu0qRJk/SjH/1IgwYN0q233qr169drz5492rZtW5P7pNYDANoSn08dv5BJkyZ5/j1o0CANHjxYvXr10rZt2zR69Ogm9Wmz2c64zhsAALQMPXv2VNeuXVVaWqrRo0fLbrfr6NGjXm3q6up0/Pjxc15aRq0HALQlfn+813eLr6QmFV8AANByffHFF/r666+VnJwsSXI4HKqoqFBRUZGnzZtvvqmGhgalpaUFa5gAAASM34M2xRcAgNalqqpKxcXFKi4uliQdPHhQxcXFKisrU1VVlR566CHt2rVLn332mbZu3apbbrlFvXv3VkZGhiSpX79+Gjt2rKZNm6Z3331X77zzjmbMmKFJkyZxx3EAQLvg86njVVVVnqPT0j+Lb0JCghISEjR37lxlZWXJbrfrwIEDevjhh89ZfJcsWaLa2lqKLwAALch7772nkSNHel7n5uZKkrKzs7V48WJ9+OGHWrFihSoqKpSSkqIxY8Zo/vz5Xqd+r1y5UjNmzNDo0aMVGhqqrKwsPf/88wHfFwAAgsHnoE3xBQCgbRsxYoQsyzrn9k2bNl2wj4SEBK1atcrksAAAaDV8DtoUXwAAAAAAzs3v12gDAAAAANCeELQBAAAAADCIoA0AAAAAgEEEbQAAAAAADCJoAwAAAABgEEEbAAAAAACDCNoAAAAAABhE0AYAAAAAwCCCNgAAAAAABhG0AQAAAAAwiKANAAAAAIBBBG0AAAAAAAwiaAMAAAAAYBBBGwAAAAAAgwjaAAAAAAAYRNAGAAAAAMAggjYAAAAAAAYRtAEAAAAAMIigDQAAAACAQQRtAAAAAAAMImgDAAAAAGAQQRsAAAAAAIMI2gAAAAAAGETQBgAAAADAIII2AAAAAAAGEbQBAAAAADCIoA0AAAAAgEE+B+0dO3ZowoQJSklJUUhIiNatW+e13bIszZkzR8nJyYqKilJ6ero++eQTrzbHjx/X5MmTFRsbq/j4eE2dOlVVVVXN2hEAAGAGtR4AgObxOWhXV1dryJAhKigoOOv2hQsX6vnnn9eSJUu0e/dudezYURkZGTp9+rSnzeTJk7Vv3z5t3rxZ69ev144dOzR9+vSm7wUAADCGWg8AQPOE+/qGzMxMZWZmnnWbZVlatGiRZs+erVtuuUWS9PLLLyspKUnr1q3TpEmT9PHHH2vjxo3as2ePhg8fLkl64YUXNG7cOD399NNKSUlpxu4AAIDmotYDANA8Rq/RPnjwoJxOp9LT0z3r4uLilJaWpsLCQklSYWGh4uPjPYVXktLT0xUaGqrdu3eftV+32y2Xy+W1AACAwKPWAwBwYUaDttPplCQlJSV5rU9KSvJsczqdSkxM9NoeHh6uhIQET5vvy8/PV1xcnGdJTU01OWwAAHCRqPUAAFxYq7jreF5eniorKz3LoUOHgj0kAABgELUeANCWGA3adrtdklReXu61vry83LPNbrfr6NGjXtvr6up0/PhxT5vvs9lsio2N9VoAAEDgUesBALgwo0G7R48estvt2rp1q2edy+XS7t275XA4JEkOh0MVFRUqKirytHnzzTfV0NCgtLQ0k8MBAACGUesBALgwn+86XlVVpdLSUs/rgwcPqri4WAkJCerWrZsefPBBPf7447riiivUo0cPPfLII0pJSdGtt94qSerXr5/Gjh2radOmacmSJaqtrdWMGTM0adIk7kIKAEALQK0HAKB5fA7a7733nkaOHOl5nZubK0nKzs7W8uXL9fDDD6u6ulrTp09XRUWFbrjhBm3cuFGRkZGe96xcuVIzZszQ6NGjFRoaqqysLD3//PMGdgcAADQXtR4AgObxOWiPGDFClmWdc3tISIjmzZunefPmnbNNQkKCVq1a5etHAwCAAKDWAwDQPK3iruMAAAAAALQWBG0AAAAAAAwiaAMAAAAAYBBBGwAAAAAAgwjaAAAAAAAYRNAGAAAAAMAggjYAAAAAAAYRtAEAAAAAMIigDQAAAACAQQRtAAAAAAAMImgDAAAAAGAQQRsAAAAAAIMI2gAAAAAAGETQBgAAAADAIII2AAAAAAAGEbQBAAAAADCIoA0AAAAAgEEEbQAAAAAADCJoAwAAAABgEEEbAAAAAACDCNoAAAAAABhE0AYAAAAAwCCCNgAAAAAABhG0AQAAAAAwiKANAAAAAIBBBG0AAAAAAAwiaAMAAAAAYBBBGwAAAAAAg4wH7ccee0whISFeS9++fT3bT58+rZycHHXp0kWdOnVSVlaWysvLTQ8DAAD4CbUeAIDz88sR7QEDBujIkSOe5e233/Zsmzlzpv7617/qtdde0/bt23X48GFNnDjRH8MAAAB+Qq0HAODcwv3SaXi47Hb7GesrKyu1dOlSrVq1SqNGjZIkLVu2TP369dOuXbt03XXX+WM4AADAMNO13u12y+12e167XC7/DBwAgADwyxHtTz75RCkpKerZs6cmT56ssrIySVJRUZFqa2uVnp7uadu3b19169ZNhYWF5+zP7XbL5XJ5LQAAIHhM1/r8/HzFxcV5ltTUVL/vAwAA/mI8aKelpWn58uXauHGjFi9erIMHD+rGG2/UiRMn5HQ6FRERofj4eK/3JCUlyel0nrNPii8AAC2HP2p9Xl6eKisrPcuhQ4f8vBcAAPiP8VPHMzMzPf8ePHiw0tLS1L17d/3pT39SVFRUk/rMy8tTbm6u57XL5SJsAwAQJP6o9TabTTabzdQQAQAIKr8/3is+Pl5XXnmlSktLZbfbVVNTo4qKCq825eXlZ73Oq5HNZlNsbKzXAgAAWgYTtR4AgLbE70G7qqpKBw4cUHJysoYNG6YOHTpo69atnu0lJSUqKyuTw+Hw91AAAIAfUOsBAPBm/NTxX/7yl5owYYK6d++uw4cP69FHH1VYWJhuv/12xcXFaerUqcrNzVVCQoJiY2P1wAMPyOFwcMdxAABaCWo9AADnZzxof/HFF7r99tv19ddf65JLLtENN9ygXbt26ZJLLpEkPfvsswoNDVVWVpbcbrcyMjL04osvmh4GAADwE2o9AADnZzxor169+rzbIyMjVVBQoIKCAtMfDQAAAoBaDwDA+fn9Gm0AAAAAANoTgjYAAAAAAAYRtAEAAAAAMIigDQAAAACAQQRtAAAAAAAMImgDAAAAAGAQQRsAAAAAAIMI2gAAAAAAGETQBgAAAADAIII2AAAAAAAGEbQBAAAAADCIoA0AAAAAgEEEbQAAAAAADCJoAwAAAABgEEEbAAAAAACDCNoAAAAAABhE0AYAAAAAwCCCNgAAAAAABhG0AQAAAAAwiKANAAAAAIBBBG0AAAAAAAwKD/YAAAAAWqtDx0/p9JeVRvrq3DFCl8ZHGekLABBcBG0AAAAfxUZ1kCQ9/UaJ9m2qMdJnVIcwbfnFzYRtAGgDCNp+ZKsolQ53MtNZdBcpPtVMXwAAoFkSO9kkSc9NGqrTXQc1u7/So1V68NVifVNdQ9AGgDaAoO0H9ZEJOmnZlPrWz6W3DHXaIVrKeZewDQBAC9L7kk5SSlywhwEAaGEI2n5Q2+lSpbuf0su39/q2ADfXsX9Ia6ZJJ78maAMAAABAC0fQ9pPD6qqPGnrotNX8oB1pVam39G3gNoVT0QEAaD5DtTnyWJVSdMxIXwCA4Atq0C4oKNBTTz0lp9OpIUOG6IUXXtC1114bzCEZ0bljhKI6hOnBV4uN9JeiY9pisyl6zTQj/UniVHQAQEC01Vqv6C7f1lJDtbm3pC02m/72WX9JVxrpk7uYA0DwBC1ov/rqq8rNzdWSJUuUlpamRYsWKSMjQyUlJUpMTAzWsIy4ND5KW35xs76pNnMX0tKjVUp/9Sk9kZGi1ITmF0xbRem3149zKjoAwI/acq1XfOq3X1if/NpId8fLPlLCxhw9v3639llHjfTJXcwBIHiCFrSfeeYZTZs2TT/5yU8kSUuWLNHrr7+ul156Sb/61a+CNSxjLo2PMlbYOneM0DcdkvSTTTWSmh/eB4RU6XWbdLTKLVN/5jjLPlHVN+WGepM6dU6SvdsVxvoDAAReW6/1ik819oV1wv/73yVjO8kd3/ybqx06fkq/2XSYu5i3NBWHjH0548HlgECLFJSgXVNTo6KiIuXl5XnWhYaGKj09XYWFhWe0d7vdcrvdnteVlZWSJJfL1eyxnKiqlsttffu/Bvrzh5hQae20q1Rx0swR8mOltXK9Y2n/+zv15VcVze7vVMVR9fzbTCWGmBmfJJ20IrTzxmcVFd/Kj3gAaPHiu16qLsndmt1PYw2xLKvZfbUFLanWS62g3tdFSPWRitv4MyPdJUpaowjt3d2g2qQUI32iecJPf6PUt36m0PrTRvttCIvUoZHPqy6ys9F+gbYkGLU+KEH72LFjqq+vV1JSktf6pKQk/f3vfz+jfX5+vubOnXvG+tRUg9/eLRhrrq9WIzfYAzi/394V7BEAgM9OnDihuDge99Qia73UDuv9PcEeAPzuhPTE7cEeBNCuXEytbxV3Hc/Ly1Nu7j9DYUNDg44fP64uXbooJCSkWX27XC6lpqbq0KFDio2Nbe5Q2wXmzHfMme+YM98xZ74zOWeWZenEiRNKSeHoYVP4s9ZL/PfhK+bLd8yZ75gz3zFnvgtWrQ9K0O7atavCwsJUXu59TW95ebnsdvsZ7W02m2w2m9e6+Ph4o2OKjY3lh9VHzJnvmDPfMWe+Y858Z2rOOJL9Ty2x1kv89+Er5st3zJnvmDPfMWe+C3StD232JzVBRESEhg0bpq1bt3rWNTQ0aOvWrXI4HMEYEgAAMIhaDwBoz4J26nhubq6ys7M1fPhwXXvttVq0aJGqq6s9dyYFAACtG7UeANBeBS1o33bbbfrqq680Z84cOZ1ODR06VBs3bjzjpin+ZrPZ9Oijj55xuhrOjTnzHXPmO+bMd8yZ75gz/2optV7i/2tfMV++Y858x5z5jjnzXbDmLMTiOSQAAAAAABgTlGu0AQAAAABoqwjaAAAAAAAYRNAGAAAAAMAggjYAAAAAAAYRtAEAAAAAMKhdBO2CggJdfvnlioyMVFpamt59993ztn/ttdfUt29fRUZGatCgQdqwYUOARtpy+DJnf/jDH3TjjTeqc+fO6ty5s9LT0y84x22Rrz9njVavXq2QkBDdeuut/h1gC+TrnFVUVCgnJ0fJycmy2Wy68sor291/n77O2aJFi9SnTx9FRUUpNTVVM2fO1OnTpwM02uDasWOHJkyYoJSUFIWEhGjdunUXfM+2bdt09dVXy2azqXfv3lq+fLnfxwlzqPe+odb7jlrvO2q976j1F69F13qrjVu9erUVERFhvfTSS9a+ffusadOmWfHx8VZ5eflZ27/zzjtWWFiYtXDhQmv//v3W7NmzrQ4dOlh79+4N8MiDx9c5u+OOO6yCggLrgw8+sD7++GPr7rvvtuLi4qwvvvgiwCMPHl/nrNHBgwetSy+91LrxxhutW265JTCDbSF8nTO3220NHz7cGjdunPX2229bBw8etLZt22YVFxcHeOTB4+ucrVy50rLZbNbKlSutgwcPWps2bbKSk5OtmTNnBnjkwbFhwwbrN7/5jbVmzRpLkrV27drztv/000+t6OhoKzc319q/f7/1wgsvWGFhYdbGjRsDM2A0C/XeN9R631HrfUet9x213jctuda3+aB97bXXWjk5OZ7X9fX1VkpKipWfn3/W9j/+8Y+t8ePHe61LS0uz7r33Xr+OsyXxdc6+r66uzoqJibFWrFjhryG2OE2Zs7q6OusHP/iB9V//9V9WdnZ2uyu+vs7Z4sWLrZ49e1o1NTWBGmKL4+uc5eTkWKNGjfJal5uba11//fV+HWdLdDHF9+GHH7YGDBjgte62226zMjIy/DgymEK99w213nfUet9R631HrW+6llbr2/Sp4zU1NSoqKlJ6erpnXWhoqNLT01VYWHjW9xQWFnq1l6SMjIxztm9rmjJn33fy5EnV1tYqISHBX8NsUZo6Z/PmzVNiYqKmTp0aiGG2KE2Zs7/85S9yOBzKyclRUlKSBg4cqCeffFL19fWBGnZQNWXOfvCDH6ioqMhzytmnn36qDRs2aNy4cQEZc2vT3n//t2bUe99Q631Hrfcdtd531Hr/C+Tv/nDjPbYgx44dU319vZKSkrzWJyUl6e9///tZ3+N0Os/a3ul0+m2cLUlT5uz7Zs2apZSUlDN+iNuqpszZ22+/raVLl6q4uDgAI2x5mjJnn376qd58801NnjxZGzZsUGlpqX7605+qtrZWjz76aCCGHVRNmbM77rhDx44d0w033CDLslRXV6f77rtPv/71rwMx5FbnXL//XS6XTp06paioqCCNDBdCvfcNtd531HrfUet9R633v0DW+jZ9RBuBt2DBAq1evVpr165VZGRksIfTIp04cUJTpkzRH/7wB3Xt2jXYw2k1GhoalJiYqN///vcaNmyYbrvtNv3mN7/RkiVLgj20Fmvbtm168skn9eKLL+r999/XmjVr9Prrr2v+/PnBHhqAVoxaf2HU+qah1vuOWt9ytekj2l27dlVYWJjKy8u91peXl8tut5/1PXa73af2bU1T5qzR008/rQULFmjLli0aPHiwP4fZovg6ZwcOHNBnn32mCRMmeNY1NDRIksLDw1VSUqJevXr5d9BB1pSfs+TkZHXo0EFhYWGedf369ZPT6VRNTY0iIiL8OuZga8qcPfLII5oyZYr+/d//XZI0aNAgVVdXa/r06frNb36j0FC+a/2uc/3+j42N5Wh2C0e99w213nfUet9R631Hrfe/QNb6Nj3zERERGjZsmLZu3epZ19DQoK1bt8rhcJz1PQ6Hw6u9JG3evPmc7duapsyZJC1cuFDz58/Xxo0bNXz48EAMtcXwdc769u2rvXv3qri42LP86Ec/0siRI1VcXKzU1NRADj8omvJzdv3116u0tNTzh4ok/eMf/1BycnKbL7xS0+bs5MmTZxTYxj9evr1nCL6rvf/+b82o976h1vuOWu87ar3vqPX+F9Df/cZvr9bCrF692rLZbNby5cut/fv3W9OnT7fi4+Mtp9NpWZZlTZkyxfrVr37laf/OO+9Y4eHh1tNPP219/PHH1qOPPtquHvdhWb7P2YIFC6yIiAjrv//7v60jR454lhMnTgRrFwLO1zn7vvZ4J1Jf56ysrMyKiYmxZsyYYZWUlFjr16+3EhMTrccffzxYuxBwvs7Zo48+asXExFh//OMfrU8//dR64403rF69elk//vGPg7ULAXXixAnrgw8+sD744ANLkvXMM89YH3zwgfX5559blmVZv/rVr6wpU6Z42jc+8uOhhx6yPv74Y6ugoIDHe7Ui1HvfUOt9R633HbXed9R637TkWt/mg7ZlWdYLL7xgdevWzYqIiLCuvfZaa9euXZ5tN998s5Wdne3V/k9/+pN15ZVXWhEREdaAAQOs119/PcAjDj5f5qx79+6WpDOWRx99NPADDyJff86+qz0WX8vyfc527txppaWlWTabzerZs6f1xBNPWHV1dQEedXD5Mme1tbXWY489ZvXq1cuKjIy0UlNTrZ/+9KfWN998E/iBB8Fbb7111t9NjXOUnZ1t3XzzzWe8Z+jQoVZERITVs2dPa9myZQEfN5qOeu8bar3vqPW+o9b7jlp/8VpyrQ+xLM4pAAAAAADAlDZ9jTYAAAAAAIFG0AYAAAAAwCCCNgAAAAAABhG0AQAAAAAwiKANAAAAAIBBBG0AAAAAAAwiaAMAAAAAYBBBGwAAAAAAgwjaAAAAAAAYRNAGAAAAAMAggjYAAAAAAAb9/zGsfDkp8YcHAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1200x900 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(2, 2, figsize=(12, 9))\n",
    "bins = np.linspace(0, 1, 20)\n",
    "\n",
    "for i, ax in enumerate(axes.ravel()):\n",
    "    x = p.posterior[\"p\"].sel({\"rating_dim\": f'{i + 1}'}).to_numpy()\n",
    "    ax.hist(x[..., 0].flatten(), bins=bins, histtype=\"step\", color=\"C0\")\n",
    "    ax.hist(x[..., 1].flatten(), bins=bins, histtype=\"step\", color=\"C1\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "dev",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.7"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
